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ABSTRACT 


A  numerical  analysis  of  the  equations  of  thermo- 
mechanical  flow  applicable  to  geodynamics  is  considered 
under  the  assumption  of  creep.  An  analytical  investigation 
of  these  equations  suggests  the  seismic  low  velocity  zone 
(L.V.Z.)  can  be  interpreted  as  a  region  of  stationary 
thermo- mechanical  coupling  between  the  tectonic  plates  at 
the  earth's  surface  and  the  mantle  convection  cells  below. 
With  the  introduction  of  a  steady  state  concentration 
gradient  of  melt  one  can  obtain  a  viscosity  distribution 
which  fits  with  the  observed  shear  wave  velocity 
distribution  in  the  L,¥„Z«  as  well  as  with  measurements  of 
seismic  attenuation.  The  introduction  of  inertial  terms  in 
the  thermal  and  mechanical  field  equations  permit  an 
explanation  of  deep  earthquakes  as  a  shear  heating 
instability  in  the  boundary  layer  at  the  upper  boundary  of 
the  subducting  plates.  The  deep  earthquakes  must  occur 
within  this  fluid  boundary  layer  in  order  to  produce  the 
observed  compressional  focal  mechanism.  This  differs  from 
suggestions  by  others  that  the  events  occur  in  the 
subducting  lithosphere.  The  predictions  of  the  theory 
developed  for  this  mechanism  of  deep  earthquakes  then 
exhibit  a  good  correlation  with  observations. 
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NOTATION 


The  following  symbols  are  used  throughout  the  thesis. 
Also  given  are  characteristic  values  and  the  references  from 
which  these  values  were  taken  (the  temperatures  are  given  in 
°C  and  all  other  values  are  in  c.g.s.  units  unless  otherwise 
specified) . 

Symbol  Value  Reference 


Density 


Heat  Capacitance 


Thermal 

Conductivity 


P  3 

c  1 . 3x1 07 

k  3  x  1  0  5 


MacKenzie 

(1370) 

MacKenzie 

(1970) 

MacKen  zie 


Activation 

Energy 


Gas  Constant 


Viscosity 

Melting 

Tem  pe nature 


E  120  Kcal/mole 

(Chapter  2) 

60  Kcal/mole 

(Chapter  3)  Griggs  6  Baker 

(1969) 

R  20  Kcal 

°/mole 

y  1022  Cathles  (1976) 


1400  (Low 

velocity  zone)  Ringwood  (1975) 


Boundary 

Temperature 


Strain  Rate 


€ 


1100 

(Lit hospher e 

Boundary)  Ringwood  (1975) 

10  i*  Carter  (1971) 


Only  c,  p  and  R  are  well  determined.  The  thermal 
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conductivity  kr  the  viscosity  p  and  their  depencence  on 
temperature  are  uncertain  (Mackenzie  1970).  The  temperature 
are  magnitudes  used  only  as  scaling  factors  and  thus  do  not 
affect  the  general  physical  outcome.  The  difference  between 
Tg  and  Tm  however  is  important  and  this  difference  is  very 
uncertain  as  is  the  strain  rate.  The  uncertainty  in  there 
values  is  a  result  of  the  fact  that  direct  measurement 


cannot  be  done  and  we  thus  must  rely  on  indirect  arguments 


CHAPTER  1 


Introduction 


1 . 1  The  Geodynamic  Problems  Under  Consideration 


The  pur 

pose  of 

this 

dissertation  is  to 

construct 

mathematical 

models  for 

two 

specific  geodynamic 

phenomena 

which  are 

generally 

assumed  to  be  problems 

n  -4r  1  »  -  **  £ 

ii  .4.  ,i>  U  jc  U 

dynamics.  The  first  is  the  coupled  thermo-mechanical 
behavior  of  the  upper  mantle.  This  can  result  from  the 
stress  induced  by  relative  motions  between  tectonic  plates 
and  the  upper  part  of  mantle  convection  cells. 
{Turcotte , D. L.  and  Oxburghr  E.R.  1972,  Torrance,  K. E.  and 
Turcotte,  D.L.  1971a, b,  Schubert  G.  et  al  1976,  Parmentier 
E.fi.  et  al  1976).  In  this  analysis  I  assume  a  constant 
stress.  This  however  is  the  same  as  the  assumption  of  a 
constant  velocity  boundary  condition  and  creep.  Creep  means 
that  the  time  derivatives  of  velocity  are  zero  and  thus  from 
the  Navier  Stokes  equation  one  can  see  that  this  is 
equivalent  to  constant  stress.  The  second  is  the  focal 
mechanism  of  deep  earthquakes  (Shaw,  H.R.  1969,  Griggs,  D.T. 
and  Baker,  D.W.  1969,  Randall  M.J,  1972). 

Of  primary  concern  in  the  discussions  of  these  two 
processes  is  the  stability  or  the  ability  of  the  very 
viscous  fluid  which  makes  up  the  earth's  mantle  to  withstand 
the  shear  stresses  imparted  to  it  without  melting.  In  the 
seismic  low  velocity  zone  the  temperature  and  pressure  is 
such  that  the  mantle  material  (essentially  iron  magnesium 
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silicates)  in  this  region  is  close  to  its  solidus.  (Ringwood 
A.E.  1975,  MacDonald  G.J.F.  and  Ness  N.F.  1961r  Wyllie  P.J, 
1971)  .  In  this  sense  the  seismic  low  velocity  zone  exists  in 
a  marginally  stable  state  (i.e.  close  to  melting)  even  when 
unstressed.  On  the  other  hand  the  subducting  plate  is  far 
from  instability  in  the  unstressed  state.  When  the  time  over 
which  this  material  can  achieve  thermal  equilibrium  is  large 
in  comparison  to  the  rate  at  which  heating  takes  place 
within  the  material,  we  may  consider  the  material  to  have  a 
low  thermal  conductivity.  P.W.  Bridgman  (1936)  found  that 
this  precise  circumstance  of  low  thermal  conductivity  and  a 
high  melting  point  relative  to  the  ambient  temperature  of 
the  experiment  caused  materials  to  exhibit  an  anomalous 
behavior  when  subjected  to  large  shear  stress  at  high 
pressures.  Here  high  pressure  means  that  the  pressure  is 
sufficiently  large  that  brittle  fracture  (breaking  of  the 
material  with  frictional  sliding  at  the  fracture  zone) 
cannot  occur.  This  is  probably  the  case  for  pressures 
greater  than  30  kb  which  is  the  pressure  at  approximately 
100  km  in  the  earth  (Griggs  D.  T.  and  Baker  D*  W.  1969) . 

A  series  of  experiments  were  performed  by  Bridgman 
(1935)  in  which  thin  disks  were  confined  between  hardened 
steel  parts.  The  disks  were  then  subjected  to  normal 
pressure  (50  kb)  and  torque  simultaneously.  Bridgman  found 
that  most  materials  deformed  smoothly  under  this  condition 
of  large  shear  stress  at  high  pressure.  This  appears  to  be 
the  case  in  the  seismic  lew  velocity  zone  due  to  the  absence 
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of  any  observations  of  sudden  stress  relief  in  that  region. 
However  in  the  case  of  materials  with  a  low  thermal 
conductivity  and  high  melting  point,  the  application  a  of 
large  shear  stress  caused  violent  snapping.  Among  the 
qualitative  effects  found  was  that  many  substances  normally 
stable  became  unstable  and  detonated,  and  conversely 
combinations  of  substances  normally  inert  to  each  other 
combined  explosivly.  Here  the  applied  force  would  suddenly 
drop  to  a  very  small  value  and  then  build  up  at  a  nearly 
steady  rate  until  another  snap  occurred.  Bridgman6 s  actual 
experiments  are  probably  not  of  much  use  for  geophysical 
purposes  since  the  strain  rates  were  much  to  large.  Thus  the 
actual  conditions  in  the  earth  should  allow  a  greater  amount 
of  fluid  behavior  This  is  important  since  SiO  showed 
rupture  only  without  flow. 

Bridgman  (1936)  and  others  (Gruntfest  I.J.  1963,  Shaw 
K.  R„  1969,  Griggs  D.T,  and  Baker  D.tf.  1969)  have  interpreted 
this  phenomenon  as  the  mechanism  of  deep  earthquakes.  They 
however  encounter  problems  in  the  application  of  this  theory 
to  deep  earthquakes.  It  is,  generally  acknowledged,  (e.g, 
Griggs  D.T.  and  Baker  D.W.  1969)  that  the  mantle  material 
has  very  little  strength  in  comparison  to  the  subducted 
slab.  Thus  in  the  above  papers  it  is  assumed  that  the  deep 
earthquakes  must  occur  within  the  subducted  slab.  However 
all  attempts  to  model  this  fluid  instability  within  the 
subducted  slab  have  failed  because  the  mantle  cannot  support 
sufficient  stress  to  cause  instability  in  the  slab.  (See 
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section  4*3  of  this  dissertion) 


The  mantle  must  yield  to  fluid  deformation  for  smaller 
values  of  stress  and  exhibit  larger  inertial  motions  ^non 
zero  time  derivatives  in  the  velocity  field)  than  the  slab 
because  of  its  lower  viscosity.  Also  the  higher  temperatures 
of  the  mantle  cause  it  to  be  more  ductile  than  the  slab 
iRanalli  G.  1974).  Thus  due  to  the  greater  strength  of  the 
slab  it  is  more  realistic  to  consider  the  motions  of  the 
slab  as  elastic  deformations  caused  by  the  viscous 
resistance  of  the  bounding  fluid.  One  then  obtains 
instability  in  the  fluid  at  the  boundary  of  the  slab 
relieving  the  elastic  deformation  of  the  slab. 


As  the  slab  subducts  the  surface  is  heated  by  thermal 
conduction  and  since  melting  of  crustal  fractions  will  occur 
crustal  material  must  mix  with  mantle  material.  Thus  the 
surrounding  fluid  in  which  the  instability  occurs  has  a 
composition  far  more  similar  to  the  crustal  material  than 
mantle  material.  This  diffusion  of  crustal  material  then 
also  causes  large  viscosity  gradients  due  to  the  changing 
chemical  composition.  This  viscosity  gradient  then  induces 
an  inertial  shear  heating  source  term  in  the  heat  balance 
eg  ua tion. 


1.2  Fluid  Dynamics  Content 


Fluid  dynamics  is  by  its 
macroscopic  phenomena  since  a 


very 

fluid 


nature  the  study  of 
is  regarded  as  a 
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continuous  medium*  One  thus  considers  even  in  the  so  called 
"infinitesimal  volume  element15  of  a  fluid,  the  resultant 
approximate  behavior  of  many  molecules  assumed  to  be 
contained  within  that  volume  (Landau  L.D.  and  Lifshitz  E.M. 
1959)  . 

The  description  of  a  classical  one-phase  fluid  is  given 
by  the  equation  of  motion  for  the  density  (conservation  of 
mass) ,  the  equations  of  motion  for  the  components  of 
velocity  (conservation  of  momentum) ,  an  equation  of  heat 
transfer  (the  equation  of  motion  for  the  entropy)  and 
supplemented  by  an  equation  of  state  relating  any  three 
thermodynamic  variables  for  a  given  amount  of  substance 
contained  in  the  system  (Putterraan  S.J.  1974), 

Of  primary  importance  in  determining  if  equilibrium 
solutions  exist  for  one  phase  fluid  flow  is  a  study  of  the 
relationship  between  the  heat  sources  and  the  heat  loss.  If 
the  heat  input  into  such  a  system  increases  at  a  sufficient 
rate  with  temperature  the  heat,  loss  from  the  system  may  not 
be  sufficient  to  allow  the  evolution  of  an  equilibrium  state 
(Landau  L.D.  and  Lifshitz  E.M.  1959).  This  condition  will  be 
called  a  "thermal  instability”  (Griggs  D.T.  and  Baker  D.R. 
1969,  Gruntfest  X.J.  1963). 

A  limited  form  of  a  theory  for  such  a  condition  was 
first  developed  for  exothermic  combustion  reactions  and 
called  the  "thermal  theory  of  explosions"  (Semenov  N.N. 
1928) .  A  quantitative  theory  was  later  worked  out  for  a 
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material  (Frank-Raraenetski  D.A.  1955) 
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sources  in  a  material  (Frank-Ramenetski  D.A.  1955).  The 
first  experimental  work  on  viscous  flow  in  this  area  was 
done  by  p.w.  Bridgman  (1935,1936,193?)  as  discussed  in  the 
previous  section.  The  conditions  under  which  stability 
(steady  state  solutions)  can  be  obtained  for  constant  shear 
stress  have  been  found  for  both  Newtonian  and  non-Newtonian 
stress  strain-rate  relations  (Gruntfest  I.J.  1963,  Griggs 
D.T.  and  Baker  D.W.  1969).  In  the  applications  to  deep 
earthquakes  in  the  above  papers  only  the  analysis  of  the 
temperature  field  and  its  behavior  relative  to  the  source 
term  driven  by  a  constant  stress  has  been  considered.  In 
Chapter  4  I  extend  this  work  by  considering  the  effect  of 
inertial  terms  in  determining  the  onset  of  instability. 


1.  .3 .  Mathematical  and  Numerical  content 

In  Chapter  2  a  variational  principle  is  used  to  study 
the  coupled  thermal  and  mechanical  field  equations  with  the 
first  order  time  derivative  terms  present.  Using  this 
variational  principle  I  consider  the  evolution  of  a  sheared 
slab  under  the  assumption  of  creep  (the  time  derivatives  in 
the  equations  of  viscous  flow  are  taken  as  zero).  The  values 
of  the  physical  constants  and  the  dimensions  are  chosen  here 
to  approximate  the  conditions  in  the  seismic  low  velocity 
zone.  An  interesting  result  occurs  in  this  analysis.  »hen 
one  assumes  a  viscosity  gradient  across  the  slab  one  obtains 
the  evolution  of  a  large  temperature  gradient  by  the 
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boundary  of  lower  viscosity.  One  thus  obtains  the  evolution 
of  very  low  viscosity  zone  near  that  boundary.  Here  it 
should  be  noted  that  one  problem  which  can  be  encountered  in 
analysing  a  highly  viscous  fluid  is  that  of  imposed 
stability  by  the  use  of  a  condition  such  as  constant  stress. 
In  using  this  condition  one  is  assuming  inertial  motions  do 
not  exist  which  is  the  very  condition  which  must  evolve  for 
instability  to  occur. 

In  Chapter  3  I  demonstrate  that  one  may  obtain  exact 
solutions  for  all  the  stress  strain-rate  relations  discussed 
in  that  chapter  if  the  physical  systems  can  be  described  by 
the  stationary  field  equations.  This  comes  about  as  a  result 
of  the  decoupling  of  the  thermo  -  mechanical  field  equations 
in  a  stationary  system. 

In  Chapter  4  I  show  that  whenever  an  intertial  term 

\ 

occurs  in  one  of  the  field  equations  the  other  must  also 
contain  inertial  terms  and  one  obtains  a  coupling  of  the 
thermo  -  mechanical  field  equations.  In  these  cases  exact 
solutions  have  not  been  found. 

In  the  actual  physical  systems  which  are  discussed  I  am 
dealing  with  non-linear  coupled  partial  differential 
equations  of  exponential  order.  Solving  such  equations 
numerically  is  in  general  very  difficult  and  in  many  cases 
impossible  iBellraan  R.  1953) .  The  method  which  I  use  for  the 
numerical  solutions  presented  here  is  the  variational 
principle  in  chapter  2.  This  method  however  has  many 
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numerical  problems.  For  example,  in  order  to  obtain  accurate 
values  for  the  integral  and  the  derivatives  a  large  number 
of  points  are  needed.  However  a  minimization  routine  is 
needed  to  find  the  extremes  of  the  functional  I(L,t) 
representing  the  time  integral  of  the  lagrangian.  Here  an 
additional  variable  is  obtained  for  each  field  equation 
represented  in  the  lagrangian  and  for  each  point  at  which 
the  minimization  is  done.  Since  the  amount  of  error  must 
increase  in  a  minimization  routine  with  the  number  of  points 
used  and  also  the  cost  increases  dramatically  with  each 
additional  variable  one  is  faced  with  severe  limitations. 
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Chapter  2 

Shear  Heating  with  the  Assumption  of  Creep 

2. 1  Physical  Assumptions 

I  define  the  state  of  motion  in  which  time  derivatives 
of  the  velocity  are  negligible  as  creep,  i.e.  the  inertial 
terms  in  the  equation  describing  the  velocity  field  are 
negligible.  In  this  chapter  I  will  consider  the  evolution  of 
a  sheared  viscous  slab  under  such  an  assumption. 

The  response  of  a  viscous  layer  to  a  constant  shear 
stress  when  the  viscosity  is  temperature  dependent  has  been 
discussed  by  Gruntfest  0563) .  A  geological  discussion  and 
extension  of  his  results  for  deep  earthquake  mechanisms  has 
been  given  by  Shaw  (1569).  A  similar  analysis  of  this 
phenomenon  has  been  given  by  Griggs  and  Baker  (1969)  using  a 

x 

more  complex  stress  strain-rate  relation  than  used  by 
Gruntfest  (1963)  and  Shaw  (1969)  .  The  Griggs  and  Baker 
(1969)  stress  strain-rate  relation  will  be  discussed  in 
chapter  3.  Here  I  will  summarize  and  extend  the  arguments  of 
Gruntfest  (1963). 

Consider  a  material  with  thermal  conductivity  k, 
specific  heat  c  and  viscosity 

p  =  PQexp [-a (T-Tfi) ] .  (2.1) 

Here  T^  represents  a  reference  temperature  (usually  the 
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ambient  temperature)  at  which  the  material  has  viscosity  yQ 
and  T  represents  the  temperature  at  any  position  in  the 
material.  The  exponential  temperature  dependence  of 
viscosity  is  reasonable  if  the  viscosity  is  controlled  by 
point -defect,  diffusion  and  T/T u-1.  Here  "a"  contains  the 
activation  energy.  Suppose  that  a  constant  shear  a  is 
applied  to  this  material.  The  rate  of  heat  production  due  to 
viscous  shear  heating  is 

oe  =  02exp[a<T-TB) }/uo  (2.2) 


provided  a  Newtonian  stress  stain-rate  relation  is  satisfied 

a  =  .  (2.2a) 

Here  l  is  the  strain  rate. 


The  heat  production  must  act  as  a  source  in  the  heat 
flow  equation.  The  temperature  is  therefore  governed  by 


pc 


dT 

dt 


2 

3  T 


dx “ 


2 

—  exp (a (T-T  ) ) 

B 


(2.3) 


where  x  is  a  co-ordinate  at  right  angles  to  the  bounding 
surfaces  of  the  layer  and  p  is  the  density. 


If  a  =0,  eg  (2,3)  is  associated  with  a  characteristic 


cooling  time 


(2.4) 


where  &  is  the  half  thickness  of  the  layer. 


%  1* 


If  k  =  0,  that  is,  if  the  system  is  thermally  isolated, 
it  is  easy  to  show  that  T  goes  to  infinity  in  a  time 


pcy 


fcM 


o 


2 

o  a 


(2 . 4a) 


A  real  physical  system  must  be  between  these  extremes. 


The  dimensionless- ratio 

2  0  2 

,  a  £  a 
A  =  , — 

yok 

of  these  two  times  represents  a  comparison  of  the  cooling 
time  to  the  heating  time  of  the  system*  Instability  occurs 
in  the  homogeneous  plane-shear  case  for  X  >.88  (Gruntfest 

I.J.  1963)  . 

Here  instability  represents  the  breakdown  of  the 
assumption  of  creep.  This  has  been  assumed  to  be  the  onset 

v 

of  an  explosive  instability  or  shear-melting  induced 
earthquake  (Gruntfest  I.J.  1963,  Shaw  H.  R.  1969,  Griggs  T. 
and  Baker  D.  W.  1969,  Nyland  E.  and  Spanos  T.J.T.  1976).  The 
validity  of  this  conclusion  is  not  clear  and  will  be 
explored  in  more  detail  in  chapters  3  and  4. 


2 . 2  A  Variational  Principle  for  the  Diffusion  Equation 

The  heat  flow  equation  (2.3)  is  a  diffusion  equation 
with  a  shear- heating  source  term.  In  this  section  I  discuss 
a  variational  principle  which  can  be  used  to  find 
approximate  numerical  solutions  of  the  diffusion  equation 
and  extend  it  to  the  full  physical  system  in  the  following 
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section.  When  considering  a  variational  principle  for  a 
dissipative  system  one  must  often  part  from  the  traditional 
methods  such  as  Hamilton* s  variational  principle.  (Morse 
P.M.  and  Feshbach  H.  1953).  There  one  is  dealing  with  a 
process  which  is  reversible  in  time.  In  dealing  with  a 
dissipative  system  without  time  symmetry,.  one  has  an 
irreversible  process,  and  thus  a  preferred  time  direction. 
This  is  the  fundamental  fact  used  in  the  construction  of 
this  variational  technique. 

Many  authors  have  studied  dissipative  systems  by 
variational  techniques.  Their  approaches  are  of  four  types: 

(1)  Some  investigators  (for  example,  Morse  P.M.  and 
Feshbach  H.  1953)  introduce  a  dual  absorption  system  and 
then  couple  the  two  systems  in  a  Hamilton  variational 
principle.  By  coupling  dissipative  and  absorptive  equations 
one  has  made  the  time  direction  once  again 
indistinguishable.  A  time  reversal  simply  switches  the 
places  of  the  absorptive  and  dissipative  fields.  There  is  no 
net  dissipation  in  the  total  system. 

(2)  Stability  analysis  can  be  formulated  as  a 
variational  principle  (Schechter,  R.S.  1967). 

(3)  A  generalization  of  Hamilton’s  principle  appears 
possible  (Djuk  D.  and  Ver jansevic ,  B.Z.  1971).  The  argument 
of  the  Lagrangian  of  the  system  is  modified  parametrically 
in  this  approach. 
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(4}  Other  variational  techniques  exist  for  irreversible 
processes  (Biot*  H.A.  1970).  It  is  possible  to  remove 
notions  such  as  a  Lagrangian  from  the  mathematics  by 
formulating  time" dependent  variational  principles  for 
dissipative  systems.  This  approach  yields  solutions  at  a 
particular  time  only.  New  calculations  must  be  made  for 
other  times. 

The  essential  feature  of  the  method  presented  here  is 
the  introduction  of  a  discrete  approximation  for  time 
derivatives.  It  is  then  assumed  that  the  state  of  the  system 
is  known  at  sometime*  t-e.  Thus,  the  variation  over  the 
field  variables  at  this  time  is  zero  and  one  may  find  the 
value  of  the  field  variables  at  the  time  t  by  demanding  that 
a  functional  L  be  stationary.  If  e  is  small  this  approach 
allows  one  to  remove  the  time  integration  from  the 
variational  principle. 

This  technique  (which  is  discussed  in  greater  detail  in 
Appendix  A)  is  illustrated  with  a  simple  example.  Consider: 

6L=6  /  —  r  (x,  t)  -<j>  (x  ,  t-e  )  ]  2+“  [  v<f>  (x , t ) ] 2  }d v  '(2.5) 

where  x  is  a  position  vector. 

The  condition  6L  =  0  will  generate  an  approximate  form 
of  the  diffusion  equation  provided  suitable  boundary  and 


initial  conditions  are  chosen. 
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After  applying  the  boundary  condition  and  (x,t-e)=0  one 

obtains : 


't'  (x,t)  }  6  <f>  (x,t)  dv 


The  term 
becomes 


in  braces  must  be  zero  and  in  the  limit  e->0  this 


a  3 1  a 


V^4)  =  0 


where  <p  =a  (T-T  )  . 


2^3  The  Coupled  Navier  Stokes  and  Heat  Flow  Equations, 


One  can  find  a  similar  variational  principle  for  the 
coupled  heat  flow  (eg  2.3)  and  Navier  Stokes  equations  for 
viscous  flow.  The  Navier  Stokes  equations  are: 


pu  -  V* (yVu)  +  VP  -  pge^  =  0 


-V 


where  u  represents  the  time  derivative  of 
field*  P  is  the  pressure  and  y  is  the 
dependent  viscosity. 


(2.6) 

the  velocity 
temperature  - 


A  variational  principle  which  will  generate  these 


equations  is: 


<5l=<5/  {——  [<J>  (x,tH  (x,  t-e)  ]  2+~  [  V  4>  (x,t)  ]  2+—  e  xp  ( —  (x,t))g2^ 

v2ae  2a  2  '  J 


^  (2.6a) 

0  [u(  x,t)-u(  x,t-eil2  +  VP.u(x,  t)  -Pge3.u  (x,t)  }dv=0 


Here  all  time  derivatives  are  given  approximately  over 
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an  interval  e.  The  coefficent  of  6 P  (x,t)  becomes  the  heat 
flow  equation  and  the  coefficient  of  <5  u  |  x  #  t )  becomes  the 
Navier  Stokes  equation.  The  variation  of  any  field  variable 
at  time  t-e  is  zero  and  the  field  variables  and  their 
spatial  derivatives  are  known  for  all  time  on  the  surface  S 
that  bounds  the  slab.  If  we  assume  one  dimensional  creep  and 
a  constant  pressure  across  the  slab  then  the  last  tv/o  terms 
in  eqn  (2.6)  vanish. 

2 . 4  Computation 

In  the  one-dimensional  analysis  I  evaluate  particular 
case  of  the  in tegral  (2 . 6a)  for  constant  stress  (i.e.  creep) 
after  interpolating  a  cubic  spline  through  the  integrand. 
The  derivatives  w.r.t.  x  are  calculated  on  this  spline.  A 
finite  number  of  values  of  <§>  and  u  which  will  minimize  the 
integrand  are  determined.  The  minimization  is  done  using 
Zangwill*s  modification  of  Powell *s  conjugate  direction 
algorithm  (Johnson  O.C.  and  Williams  L.L.  1973).  This 
routine  was  called  from  the  International  Mathematical  and 
Statistical  Library  edition  4  ?  available  on  the  University 
of  Alberta  computer.  The  accuracy  of  this  process  was  tested 
by  the  addition  of  points  until  the  structure  of  the 
velocity  and  temperature  profiles  were  smooth  across  the 
slab.  For  the  perposes  of  this  analysis  only  approximate 
values  of  the  temperature  and  velocity  are  necessary. 

Here  6  L=0  does  not  necessarily  imply  a  unique 
solution.  The  functional  L  contains  a  term  of  more  than 
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The  uniqueness  of  the  solutions  here  may  be 
tested  numerically.  This  is  done  by  making  small 
perturbations  in  the  initial  guesses  and  changes  in  the 
convergence  criteria  to  see  if  the  same  solution  will 
result. 

2^5  Some  Numerical  Examples 

In  figure  (2.1)  the  values  of  the  lagrangian  specified 
in  eqn.  (2.5)  are  shown  as  a  function  of  temperature  for 
different  values  of  the  dimensionless  constant  lambda.  Here 
a  constant  stress  is  assumed  and  the  time  step  is  kt/p  c^2  =  1 . 
In  figure  (2,2)  the  solutions  of  the  differential  equation 
(2.3)  are  illustrated  on  contour  plots  of  the  lagrangian 
specified  in  eqn.  (2.5)  with  time  steps  of  kt/pe£2=. 25, . 5, . 75 
and  1.  Here  temperature  is  plotted  on  the  vertical  axis  and 
lambda^  on  the  horizontal  axis. 

The  values  of  the  specific  heat  c,  thermal  conductivity 
k  and  the  activation  energy  for  point  -  defect  diffusion 
vary  widely  for  geologic  materials.  In  the  calculations  in 
this  section  the  values  c=.25  cal/g/C0;  k=.03  c a 1/C ° /era /sec ; 
a=.01  /C°  are  used. 

Although  these  numbers  may  not  fit  any  particular  rock 
exactly  they  do  not  vary  by  more  than  a  few  orders  of 
magnitude  for  most  geologic  materials.  The  density  p  can  be 
taken  as  3  g/cra3  within  10%  for  most  geologic  materials*  at 
”low"  pressures,  but  the  viscosity  can  vary  widely. 


-  } 


17 


If  we  use  a  Newtonian  stress  strain-rate  relation  with 
a  viscosity  of  1022  poise  and  a  velocity  gradient  of 
10“H/5ecf  the  assumption  that  thermal  instability  occurs 
leads  to  a  half  thickness  £  in  the  constant  A  of 
approximately  6x106  era.  This  value  must  be  slightly  small 
since  the  actual  value  of  A  must  be  reduced  due  to  shear 
heating.  This  crude  approximation  thus  yields  a  value  of  £ 
which  is  fairly  close  to  the  thickness  of  the  seismic  lew 
velocity  zone.  It  is  thus  possible  that  the  seismic  low 
velocity  zone  is  a  marginally  stable  material.  This  would 
seem  to  fit  well  with  the  presence  of  melt.  However  if  a 
"marginally  stable”  material  is  able  to  exist  in  the  seismic 
low  velocity  zone  without  experiencing  any  instabilities 
(observed  sudden  relief  of  stress)  one  roust  then  question 
the  validity  of  assuming  the  breakdown  of  creep  to  be  a 
representation  of  coupled  thermo  -  mechanical  instability. 
This  problem  will  be  discussed  in  both  chapters  3  and  4« 

Consider  now  some  numerical  results  assuming  creep. 
First  assume  we  have  a  slab  of  viscosity  1022  poises 

throughout.  The  slab  has  a  half  thickness  of  60  km  and  the 

gradient  of  the  velocity  field  is  given  by  o /U 
throughout  the  slab.  The  stress  within  the  slab  is  assumed 
to  be  100  bars  (the  values  of  stress  are  chosen  to  give  a 
velocity  difference  of  10"7  cm/sec  accross  the  slab) 

everywhere  and  in  all  cases  the  boundary  temperatures  are 

fixed.  I  am  therefore  considering  here  whether  the  decrease 
in  viscosity  in  the  seismic  low  velocity  zone  can  be  solely 
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attributed  to  shear  heating.  The  evolution  of  the 

temperature  and  displacement  is  shown  in  figures  (2.3  and 

2.4)  for  a  time  step  of  1014  seconds.  It  should  be  noted 
that  the  viscosity  does  not  decrease  by  the  almost  2  orders 
of  magnitude  necessary  to  be  in  agreement  with  the  observed 
viscosities  in  this  region  (see  section  3. 1  of  this 
dissertation) . 

Increasing  pressure  with  depth  and  the  rise  of  melt 
towards  the  surface  results  in  a  viscosity  which  increases 
with  depth.  The  evolution  of  a  slab  under  this  condition 
from  a  constant  temperature  throughout  may  be  illustrated  by 
y0  ->P0  exPi2(x+1)).  The  results  of  this  evolution  of 

temperature  are  then  given  in  figure  (2,5),  The 

corresponding  evolution  in  velocity  is  given  in  figure 

(2.6) ,  The  earth  also  has  a  temperature  which  increases  with 
depth  x due  to  cooling  at  the  surface.  The  evolution  of  the 
same  slab  as  in  figure  (2.5)  only  from  a  linear  temperature 
gradient  of  3C0°K  across  the  slab  is  illustrated  in  figure 

(2.7) .  The  corresponding  evolution  in  velocity  for  this  slab 
is  given  in  figure  (2.8).  Another  example  of  interest  to  the 
analysis  in  chapter  4  is  that  of  a  temperature  increase 
towards  the  boundary  of  low  viscosity.  The  evolution  of 
temperature  and  velocity  in  this  case  is  illustrated  in 
figures  (2.9)  and  (2.10)  respectively.  Figure  (2.11)  illus¬ 
trates  a  proposed  viscosity  depth  profile  in  the  seismic  low 
velocity  zone. 
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of  the  lagrangian  in  sec  (2* 2)  are  shown 
function  of  temperature  for  different 
of  the  dimensionless  constant  lambda. 
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Figure  2.2:  Solutions  of  the  variational  principle  in 

section  2.2  are  illustrated  for  4  different 
sized  time  steps,  .25,  .5,  .75  and  1..  The 

initial  condition  is  a  constant  temperature 
across  the  slab. 
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Figure  2.3:  The  evolution  of  the  temperature  field  is  shown 

from  an  initially  constant  value  across  the 
slab.  The  fluid  is  homogeneous  and  10  time 
steps  of  1 0 1 3  sec  (i.e.  episilon  =.111)  are 
used  in  the  calculation.  The  temperature  change 
is  given  in  10-~2  x  °C.  Here  the  stress  is  50 
'bars , 
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Figure  2.4:  The  velocity  evolution  is  given  corresponding 

to  the  temperature  evolution  of  the  homogeneo 
slab  in  figure  (2.3).  The  values  on  t 
contours  are  the  velocities  in  cm/sec  relative 
to  the  centre  of  the  slab. 
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Figure  2. 


>:  The  temperature  evolution  with 

gradient  P0=  U0exp  3  (Z-t- 1 ) /2)  is 
Here 
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a  viscosity 
illustrated  * 


the  stress  is  10  bars  and  the  time  steps 
10s3  sec.  The  initial  temperature 


distribution  is  constant  across  the  slab. 
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Figure  2*6: 


The  velocity  evolution  correspon 
temperature  distribution  in  figure 
inhomogeneous  slab  is  given.  The 
low  viscosity  is  on  the  right. 
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Figure  2.7:  The  earth  has  a  known  temperature  gradient. 

This  figure  illustrates  the  temperature 
distribution  for  the  slab  in  figures  (2.5)  and 
(2.6)  with  the  temperature  of  the  boundary  of 
high  viscosity  held  at  300°C  greater  than  the 
boundary  of  low  viscosity.  Here  the  stress  is 
10  bars. 
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Figure  2.8: 


The  velocity  evolution  corresponding  to  the 
slab  described  in  figure  (2® 7)  is  given.  Here 
one  obtains  a  fairly  uniform  velocity  gradient 
across  the  slab. 
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Figure  2.9:  A  temperature  increase  of  60°C  at  the  boundary 

of  low  velocity  over  the  boundary  of  high 
viscosity  is  considered.  This  figure 

illustrates  that  the  shearing  of  the  slab  then 
causes  the  maximum  temperature  to  move  awa 
from  the  boundary  of  low  viscosity  towards  th 
interior.  Here  the  time  steps  in  1013  sec  an 
the  stress  is  1  kb  and  the  slab  width  is  5  km. 
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Figure  2.10:  The  velocity  evolution  corresponding  to  the 

temperature  distribution  in  figure  { 
given.  Note  that  the  maximum  velocity  g 
occur  in  the  region  of  maximum  temperature. 
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Figure  2.11: 
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Chapter  3 


Exact  Solutions  for  Constant  Stress  and  Their 
Application  to  the  Seismic  Low  Velocity  Zone. 

3. 1  Physical  Discussion 

It  is  generally  acknowledged  (Schubert  G.  et  al,  1976) 
that  many  geophysical  characteristics  like  topography, 
lithosphere  thickness  and  heatflow  on  ocean  bottoms  are 
essentially  a  function  of  the  age  of  the  ocean  bottom. 
Within  limits  these  features  can  be  predicted  as  a 
consequence  of  relatively  simple  models  of  heat  flow  and 
dynamics  in  the  earth  (Turcotte  D. L.  and  Oxburgh  E. R.f  1967* 
McKenzie  DcP.r  1967*  Parker  R. L.  and  Oldenburg  D.W.,  1973 
for  example) .  With  rare  exceptions  (Schubert  G.  and  Turcotte 
D«L.f  1972  Schubert  G.  et  al,  1976)  none  of  these  models 
include  the  thermo- mechanical  coupling  that  must  exist  if 
shear-heating  contributes  to  temperature,  and  viscosity  is 
temperature  dependent. 

All  plausible  relations  between  stress  and  strain  rate 
in  rocks  under  mantle  conditions  involve  the  dissipation  of 
heat  in  some  form.  These  mechanisms  of  deformation  are 
controlled  by  the  migration  of  crystal  defects  and  as  such 
require  the  generation  of  heat.  This  heating  modifies  the 
physical  properties  of  the  material  and  as  a  result  can 
alter  its  mechanical  behaviour.  Under  most  geologic 
conditions  such  thermal  feedback  appears  to  be  a  fairly 
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small  contribution  to  geodynamics  but,  with  sufficiently 
effective  heat  generation  mechanisms,.  instabilities  can 
arise.  The  conditions  under  which  the  feedback  mechanism 
becomes  important  to  the  geodynamic  models  being  considered 
is  of  fundamental  importance  both  in  this  chapter  and 
chapter  4. 

Here  lithosphere  formation  will  not  be  dealt  with 
directly.  Instead  in  this  chapter  I  attempt  to  clarify  the 
part  shear  heating  plays  in  the  existence  of  the  seismic  low 
velocity  zone.  The  results  of  this  analysis  are  dependent  on 
the  many  assumptions  made  about  the  physical  parameters  and 
constants  involved*  For  example  it  is  assumed  here  that  the 
seismic  low  velocity  zone  behaves  essentially  as  a  viscous 
fluid  over  the  time  scales  of  the  motion  imparted  to  it.  The 
mantle  material  below  is  then  able  to  exhibit  solid  behavior 
relative  to  the  relaxation  time  for  the  seismic  low  velocity 
zone.  Yet  this  material  beneath  the  seismic  low  velocity 
zone  is  able  to  behave  as  a  fluid  itself  over  longer  periods 
of  time.  In  reality  this  approximation  represents  a  fluid 
whose  viscosity  increases  with  depth.  Highly  viscous  fluids 
do  have  the  ability  to  act  as  solids  or  fluids  for  different 
rates  of  application  of  stress  (Kolsky  H.  1963) .  It  appears 
however  that  this  is  not  a  factor  in  the  geodynamic 
processes  being  considered.  Cathles  (1975)  claims  that  there 
is  no  evidence  for  a  stress  rate  dependent  viscosity  in 
glacial  rebound  studies. 

Certainly  the  most  controversial  physical  parameter  is 
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in  fact  the  value  of  viscosity.  All  of  the  present  knowledge 
of  the  viscosity  of  the  mantle  comes  from  glacial  rebound 
studies.  A  fairly  complete  review  of  this  subject  and  list 
of  references  is  given  by  cathles  (1975).  It  should  be  noted 
however  that  Cathles  assumes  a  Newtonian  mantle.  This  is  an 
assumption  which  is  not  generally  accepted  (e.g.  Carter  n.L. 
1976r  Froidevaux  C«  and  Schubert  G.  1975) .  There  are 
differences  in  the  rates  of  change  of  viscosity  with  depth 
for  Newtonian  and  non- Newtonian  flow  lav/s.  and  the  non- 
Newtonian  (power)  flow  laws  give  rise  to  a  well  defined  low 
viscosity  channel  in  contrast  to  the  linear  law  (Carter  N.L. 
1976)  .  Cathles  (1975)  however  claims  that  the  load  cycle 
behavior  of  earth  models  that  have  high  viscosity  lower 
mantles*  non-Newtonian  mantle  viscosity  or  substantial  non- 
adiabatic  density  gradients  are  similar  to  the  load  cycle 
behavior  of  layered  Newtonian  viscosity  models.  Thus  for  the 

N. 

purposes  of  this  analysis  I  will  assume  the  values  of 
viscosity  obtained  fcv  Cathles. 

The  areas  of  uplift  which  have  been  studied  are 
Fennoscandia *  the  Canadian  Shield,  Lake  Bonneville, 
Greenland  and  the  Arctic.  The  area  over  which  uplift  occurs 
is  used  along  with  the  rate  and  horizontal  variations  in  the 
rebound  to  obtain  models  for  the  values  of  the  viscosity  to 
different  depths  under  each  region.  The  general  values 
obtained  by  cathles  are,  a  region  of  relative  low  viscosity 
of  4x1 02  0  poise  from  the  lit hosphere-a sthenosphere  boundary 
to  a  depth  of  approximately  75  km  below  the  plates  and  a 
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viscosity  of  approximately  1022  poise  throughout  the  rest  of 
the  mantle®  Although  these  values  are  subject  to  some 
dispute  (eg.  MacDonald  G.J.H.  1963)  they  probably  cannot  be 
very  much  in  error  above  the  650  km  discontinuity.  If  the 
viscosity  were  much  less  than  this  glacial  rebound  would 
have  to  occur  at  a  faster  rate  than  observed.  For  much 
larger  values  of  viscosity  one  would  be  required  to  alter 
the  equations  used  by  Cathles  governing  the  dynamics  in 
order  to  obtain  a  reasonable  rate  of  rebound.  Such  an 
analysis  certainly  does  not  seem  to  be  justified  since 
similar  values  were  obtained  for  all  areas  which  have  been 
studied  in  spite  of  the  different  areal  extents  and  load 
releases  involved. 

The  values  of  the  other  physical  constants  used  in  the 
equations  of  heat  flow  and  fluid  flow,  specific  heat  c , 
density  p  and  thermal  conductivity  k,  seem  to  be  fairly 
consistent  throughout  the  literature  on  heat  flow  in  the 
mantle,  (eg  McKenzie  D.P.  1970,  Toksoz  N.M.  et  al  1971).  The 
values  of  these  constants  are  generally  assumed  to  he 
c~ 1 07er gs/°C/gm ,  p=3gra/cia2#  k=4x 1 Osergs/cm/°C/sec. 

3® 2  A  Simple  Newtonian  Stress  Strain- Fate  Relation 

Consider  now  the  solutions  which  "may”  evolve  for  an 
infinite  slab  of  half  thickness  £  subjected  to  a  shear 
stress  a  .  The  rate  of  heat  production  due  to  viscous  shear 
heating  at  any  point  in  the  slab  is  oe  where  £  is  the 
strain  rate  and  a  is  the  stress  at  that  point.  The 
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boundaries  of  the  slab  are  held  at  temperature  T  «  The  slab 

B 

has  thermal  conductivity  k,  specific  heat  c  and  viscosity 
y  =  PQexp  (-a  (T-Tg  ))  .  Here  Pq  is  the  viscosity  at 

temperature  Tg  and  a  is  a  material  constant.  Thus  the  rate 
of  heat  production  is  given  by  eg.  (2.2)  if  a  newtonian 
stress  strain  rate  relation  is  satisfied  and  the  principle 
of  heat  balance  is  once  again  given  by  eg.  (2.3).  Momentum 
balance  eg  (2.5)  is  conveniently  derived  from  the  condition 


3£  =  d  (pv) 
3x  dt 


(3.1) 


where  v  is  the  velocity  in  the  plane  of  the  slab.  The 
complete  solution  of  a  sheared  slab  requires  simultaneous 
solutions  of  equations  (2.2) ,  (2.3)  and  (3.1). 

The  assumption  of  Gruntfest  (1963) ,  Griggs  and  Baker 
(1969) t  Shaw  (1969)  ,  Schubert  et  al  (1976)  and  the  previous 
chapter  that  (  9<j  /g  x)  =0.  But  when  inertial  term  exist  in 
the  temperature  field  this  is  not  rigorously  valid.  In  the 
analysis  of  the  inertial  field  equations  given  in  Chapter  4 
it  is  shown  that  the  stress  is  uniform  throughout  the  slab 
if  the  inertial  terms  (derivatives  w.r.t.  time)  are  zero  in 
both  the  heat-flow  and  Navier  Stokes  equations.  I  do  not 
mean  to  imply  that  this  assumption  of  creep  is  not  a  useful 
approximation  in  many  cases.  However  I  do  wish  to  point  out 
that  in  analyzing  a  system  changing  in  time  the  creep 
assumption  can  ''impose"1  stability  on  highly  unstable 
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physical  situations  and  it  must  always  breakdown  in  the 
"neighbourhood  of  instability". 

The  fashion  in  which  stability  can  be  "imposed"  is  as 
follows.  It  will  be  shown  in  chapter  4  that  a  thermal 
instability  is  a  localized  phenomenon  within  a  fluid. 
However  the  instability  criterion  derived  in  chapter  2 
depends  only  on  the  slab  as  a  whole.  Keeping  these 
statements  in  mind  it  can  be  seen  from  eq  (3.1)  that  any 
breakdown  in  the  "assumption  of  creep"  can  induce  a  large 
local  velocity  gradient  without  a  substantial  reduction  in 
stress.  This  may  then  produce  a  local  weakness  within  the 
slab  from  which  it  is  not  able  to  recover  before  the  onset 
of  melt.  If  one  were  to  then  take  a  slab  which  behaved  in 
this  fashion  and  apply  the  creep  assumptions  to  it,  the 
velocity  field  would  be  averaged  out  over  the  slab  as  a 
whole  and  the  slab  could  appear  very  stable.  In  the 
neighbourhood  of  instability  the  creep  assumptions  must 
breakdown  by  definition  since  a  sufficient  reduction  in 
viscosity  must  result  in  the  relief  of  stress  by 
accelerations.  The  analysis  of  these  phenomena  will  be  the 
subject  of  chapter  4. 

The  rest  of  this  chapter  is  concerned  strictly  with 
analytic  solutions  to  the  stationary  field  equations  and  the 
applications  discussed  in  section  (3.1)  to  the  seismic  low 
velocity  zone. 


In  general  one  has  the  condition 
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4 


w* 


cf2  (x,y) 


(3.2) 


exp (a (T-Tb) ) 


for  the  heat  balance  in  an  infinite  fluid  layer  in  which 
flow  is  not  restricted  to  be  laminar.  Here  x  is  the  co¬ 
ordinate  across  the  slab  and  y  is  a  co-ordinate  in  the  plane 
of  the  slab  and  the  direction  of  the  applied  shear.  If  there 
is  no  stress  in  the  x  direction  then  vrx,y)=0.  Further,  if 
we  assume  that  the  temperature  is  time  independent  at  any 
point  in  the  slab  then  the  velocity  must  also  be  time 
independent  and  must  be  a  function  of  y  only.  Equation  (3,2) 
then  becomes 


exp (a (T-Tb) ) 


If  the  stress  varies  in  the  horizontal  direction  one  obtains 
a  stress  gradient  and  thus  from  eqn.  (3.1)  we  have  inertial 
motions.  However  if  the  flow  is  incompressible  such  a 
condition  of  a  constant  stress  across  the  slab  and  a 
variable  stress  in  the  plane  parallel  to  the  slab  is 
impossible.  Also  if  shear  stress  causes  heating  and  in  some 
cases  melting  it  must  diffuse  at  a  finite  rate  which  may  be 
zero  across  100fa  melt.  Therefore  3a  /3  y#0  implies 
3a  /3  x*0.  I  therefore  consider  only  the  fully  stationary 
solution,  a  not  a  function  of  x,  y,  or  t.  The  inertial  term 


in  the  thermal  equation  must  therefore  be  zero  and  this 
impl  ies  3T/3y  =  0  for  \ri-  0. 
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The  differential  equation  for 


the  stable  temperature 
distribution  which  may  evolve  is  qiven  by 


exp  [a  (T-TJ  ]/p 

O 


(3.3) 


with  boundary  conditions  T  =  TB  for  x=^±  £  where  2  £  is  the 
thickness  of  the  slab. 

The  solution  of  equations  (2»2f  3 .  1 ,  3.3)  is  similar  to 
one  in  Landau  and  Lifshitz  (1959,  p«  190).  In  order  to 
reduce  the  calculation  to  dimensionless  quantities,  as  in 
chapter  2,  let 


=  a  (T-T^)  ,  and  £  =  x/£ 


d  4>  ,  2  2 

Then  - o  +  AexP(<M  =  0  where  >  = 

(3.4) 

Equation  (3.3)  has  the  solution 


(3.5) 


where  <})  is  the  maximum  temperature  (which  must  occur  at 
the  centre  of  the  slab) .  I  obtain  ^  from  the  boundary 
condition  <S>  =0  at  £  =1  and  observe  that  solutions  of  this 
equation  for  <S>  exist  only  if  \  <.88  as  illustrated  in 
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figure  (3.1).  Figure  (3.2)  shows  graphical  presentations  of 
this  solution. 

From  equations  (3.1)  and  (3.3)  one  finds  the 
nonvanishing  component  of  velocity  is  given  by 

_  _  k  dT  ka  3tj> 

a  dx  o£ 

=  sgn(^)~  exp  (<f>c/2)  tanh  f  exp  (<J>c/2  |  £| ) 


Figure  (3.3)  show  graphical  representations  of  the  velocity. 
Note  that  the  critical  parameter  is  the  value  of  A  .  For  the 
particular  case  solved  here  the  condition  that  stable 
solutions  may  exist  is  ^<*88.  It  is  however  interesting  to 
consider  the  case  where  melting  occurs  at  a  finite 
temperature  T^.  Then  the  time  to  melting^  ,  as  defined  in 
chapter  2  becomes 

t  =  _ 1 _ 

M  2  (1-exp (-a (T  -T  ) ) ) 


a  a 


M  B' 


For  the  low  velocity  zone  a(TM~TB)-1  hence  tc/tM  -.6  A  ,  but 

from  figure  (3,1)  it  can  be  seen  that  the  value  of  A 

changes  very  little  to  .85  for  this  value  of  <p  .  It  should 

c 

be  noted  that  in  this  case  t^/t^A  since  t  is  the  time  for 
the  material  to  melt  and  A  is  defined  under  the  assumption 
that  the  material  could  be  heated  to  an  infinite  temperature 
without  melting. 
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Figure  3.1: 


The  central  temperature 
dimensionless  constant 
that  the  maximum  value 
solutions  exist  is  .88 
temperature  of  1.2. 


as  a  function  of  the 
lambda  is  plotted.  Note 
of  lambda  for  which 
which  corresponds  to  a 
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Figure  3. 


h  The  values  of  the 
the  middle  x=0,0 
class  of 
stationary 
represents 
dimensionless 
are  given, 
temperature  $> 


temperature  are  given  between 
and  the  boundary  x=  1 . 0  of  the 
homogeneous  slabs  in  which  a 
solutions  exist.  The  vertical  axis 
different  values  of  the 

constant  lambda.  The  temperatures 
in  terms  of  the  dimentionless 
lie.  10-2  x  oC) e 
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Figure  3.1 


:  The  values  of  the  velocity  are  given 

same  slabs  as  in  figure  (3.2),  The 
the  velocity  are  on  the  contour 
cm/sec* 
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i«.3  Solutions  for  more  Realistic  Stress  St  rain -rate 
Relations 

A  further  complication  in  the  geophysical  application 
of  these  equations  is  the  necessity  of  using  "realistic” 
stress  strain  rate  relations.  This  is  a  subject  which  has 
been  treated  very  extensively  by  Ranalli  ( 19 74 )  and  mere 
recently  in  a  review  paper  by  Carter  (1976)  and  thus  the 
justification  of  the  stress  strain-rate  relation  will  not  be 
considered  here. 

Steady  state  creep  rates  for  metals,  ceramics  and  rocks 
are  related  to  temperature  and  stress  by  a  Weertman  -  type 
equation  of  the  form  (Carter  N.L.  1976) 

e  =  A  exp(-(E+PV)/RT)f  (cr)  (3.6) 

\ 

Here  A  is  a  slightly  temperature  sensitive  material 
constant,  E  creep  activation  energy,  P  is  the  pressure,  V  is 
the  activation  volume,  R  is  the  gas  constant,  T  it  the 
temperature  and  f(o  )  is  the  stress  function.  In  this 
section  I  will  be  concerned  mainly  with  the  special  case  of 
eqn  (3. 6) , 
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^  ~  exp  j 
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RT  / 
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) 

(3.7) 

which  holds  for  most  geologic  material  in  the  neighbourhood 
of  instability  (Griggs  D.T.  and  Baker  D.W.  1969). 

Writing  this  relationship  in  a  form 


which 


includes 
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,s viscosit y"  one  has 


n  n 


o  =  Vie, 


o 


o  (  E 

y  =  —  exp  j 


E 


o 


l  nRT  nRT 

\ 


B 


Where  n  is  an  experimentally  determined  constant.  Observe 
that  for  small  temperature  changes  (T  ^TB)  we  can  write 


exp 


E  f  1  1  \ 


nR  \  Tb  T 


~  exp  1 


E 


V  nRT 


T  {VT> 


(3.8) 


B 


so  that  the  constant  ’^a11  used  previously  becomes 


E 


nRT 


B 


and  y0  becomes  oQ/  e0. 

The  approximation  given  by  equation  (3.8)  relies  on  the 
condition  that  at  least  in  the  seismic  low  velocity  zone  the 
actual  temperature  is  close  to  the  melting  temperature.  This 
is  the  essence  of  a  stability  criterion  developed  by 
Robertson  (1948).  If  oQ  is  taken  as  o  the  constant  stress 
in  the  slab,  the  steady  state  heat  flow  becomes 


Ek  d'  $ 
R  dx2 


k exp  [k -  * 1 0 


where  <$>  =  (R/E)  T.  Observe  that  the  source  term  depends  on 
temperature,  on  the  constant  applied  stress  a  and  yQ  which 
is  also  constant.  This  equation  is  the  analogue  of  equation 
(3.4).  Introducing  the  dimensionless  length  £  =x/ £ 


d  cp 

<H2 


R£  i  1  12 

-  exp  —  ~  -r  o 
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which  becomes 


-j 2  , 

+  A  '  exp  ( — 1/4) )  =  0 

d£ 


(3.9) 


where 


A ' 


Ekp 


a2  exp  ( 1/ 4>B )  =A  ^B2exp  (1/4>b)) 


(3.10) 


Here  A  was  defined  in  (3.4)  and  a  was  derived  above  for  the 
assumption  T  - TB»  Multiplying  both  sides  of  (3.9)  by  d$  /d£ 
and  integrating  one  obtains 


■l/4> 


dd) 

d£ 


=  ±/2Tr 


"(j)e 


"Ei 


$  J  j 


4> 


J 


4) 


c 


where  Ei(x)  = 


x 


—  CO 


—  dt  r  is  the  non-dimensional 
tt  c 


temperature  at  the  centre  of  the  slab. 
Define  a  function  A  by 

A  ( (■> )  =  f  ~(j)e”1//^)-Ei  (-!/({))  1 


'B 


Then 


d$ 


$  [A((|,)-A  (cf>c)  ] 


X  =  ±/2X 


(3.11) 


where  <»>  the  value  of  4)  at  £=1, 
13 

determines  $  and 
c 


d<{> 


$ 


IA(4>) -A(4>c)  ] 


■y-  =  ±£/2Tt  determines  <H£) 


(3.11a) 


c 


58 


A  plot  of  i>  (X  pxj  appears  in  Figure  (3,4).  The  velocity 
is  given  by 


v 


_M  (2/A(A(cMO  -A(<f> J))"2 

Eai  8£  PQ  c 


where  the  sign  on  the  square  root  is  chosen  to  correspond  to 
positive  or  negative  E,  .  Compare  this  with  equation  (3.5)  . 

If  $  =<f>c  -  e  where  e  /<pc  «1  then 
(A  (<$>  )  -  A  ($>c  )  )  *  '2  =  (  e  exp  )  }  *  /2 

and  thus 


We  may  then  use  the  following  approximation  in  equation 


(3,5) 


tanh 


r 


a 


exp 


(  ♦c-*bV2+] 
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J  J 


exp 


24> 


2 

B  J 


exp 


1 


to  obtain 


v  ~ 


lo 

V> 


K  exp 


C 


The  introduction  of  a  melting  temperature  is  now  natural.  If 
we  argue  that  $  where  <$>M  is  the  melting  temperature 
it  is  easy  to  use  (3.11)  to  evaluate  X  at  the  condition 
@c  =$>M  .  Further  it  is  worth  noting  that  stable  solutions 
exist  for  all  values  of  X  .  Thermal  runaways  do  not  occur  in 
the  same  way  as  in  the  previous  case.  Mechanical  instability 
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in  the  material  described  by  (3.7)  comes  about  only  by 
melting . 

Now  consider  the  problem  of  determining  a  stability 
criterion  for  a  slab  heated  from  within  by  viscous  shear 
heating  with  the  stress  strain  rate  relation  (3.7).  If  we 
assume  zero  thermal  conductivity  and  uniform  shear 
throughout  the  slab  then  the  centre  will  be  heated  to  the 
melting  temperature  in  a  time 
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no  heat  sources  the  slab  has  a  characteristic  cooling 


By  analogy  with  Gruntfest  (1963),  let 

G  =  !e  =  *'/f<vV 
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where  A*  is  defined  by  equation  (3.11)  or  by  equation 
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If  tM  is  sufficiently  large  compared  with  t  it  is 
clear  that  stability  exists*  The  instability  occurs  when  the 
central  temperature  in  the  slab  reaches  the  melting 
temperature.  The  corresponding  value  of  X  is  the  critical 
value  and  can  be  calculated  from  equation  (3.11a).  It  is 
easy  to  show  that 
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is  the  limiting  value  of  G  for  stability.  When  <f>_  =  <f>  -e 


where  <<1 
'"A 


q  =  2. 

Another  stress  strain-rate  relation  of  considerable 
interest  for  steady  state  creep  in  rocks  is  (Weertman  and 
Weertman,  1375) 

e  =  Cnonexp(-E/RT)/T 

Here  the  value  of  A  from  eqn  (3 . 6)  has  a  1/T  dependence. 

The  values  of  E,  C  and  n  are  subject  to  considerable  debate, 
but  the  analytic  properties  of  flows  which  couple  to  thermal 
behaviour  are  interesting. 

Solving  for  the  heat  source  term  we  get  for  constant 
stress 

a£  =  cn  n+1  exp (-E/RT) 
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Setting 


we  get 


3%  ,  y  exp  (-1/4)) 
352  + 


=  0 


where 
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can  be  solved  numerically  for  $  .  A  stability  criterion  can 
be  derived  exactly  as  was  done  for  the  previous  relation. 

The  preceding  arguments  involving  non-linear  rheology 
are  hard  to  interpret  in  an  immediately  useful  geophysical 
way.  The  success  of  a  simple  Newtonian  stress  strain  rate 
relation  in  geodynamics  suggests  we  should  at  least  look  at 
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the  way  t  her  Biomechanical  coupling  can  modify  apparent 
viscosity*  This  theoretical  development  will  now  be 
completed  with  some  connections  to  Newtonian  processes  in 
which  thermo-mechanical  coupling  is  ignored* 


Traditionally  experiments  on  the  shear  of  layers  are 
interpreted  in  terms  of  such  Newtonian  viscosity  laws.  If 
shear -heating  modifies  the  viscosity,  the  observed 
experimental  results  are  obviously  not  a  good  measure  of 
actual  viscosity.  Calculation  of  "apparent  viscosity"  is 
fairly  easy  with  the  results  derived  above.  Here  apparent 
viscosity  is  the  ratio  of  applied  stress  to  apparent 
velocity  gradient.  V?e  can  estimate  this  quantity  using 
equation  (3.  7)  which  yields 
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(3,12) 


It  is  generally  accepted  (e.g.  Ringwood  A.E.  1975)  that  the 
temperature  of  the  seismic  low  velocity  zone  is  near  its 
solidus.  Although  there  is  no  veil  defined  melting  point  for 
such  material,  it  is  still  useful  to  introduce  this 
approximation  since  it  behaves  in  the  presence  of  partial 
melt  as  a  material  very  near  instability. 
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and  from  (3.11a) 
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Then  we  have 
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which  can  be  evaluated  as 
Stegun,  1965)  .  Observe 
fashion  on  the  stress.  In 
truncating  we  get 


exp  (\p2 )  dip  (3.13) 
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a  Dawson  integral  (Abramowitz  and 
that  y  depends  in  a  non  linear 
fact,  by  expanding  exp(i|;2)  and 
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The  stress  strain  law  used  to  derive  (3.13)  is  more 
realistic  than  the  assumption  that  viscosity  decays 
exponentially  with  temperature.  Nevertheless  another 
apparent  viscosity  can  be  derived  from  the  high  temperature 
limit  of  equation  (3.7). 


o 


y  v 
o 

l 


exp  (-<f>c/2) 


from  which  it  is  easy  to  see 
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y  =  p 
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2  exP(“^c/2) 


In  this  case  y  varies  as  o  2  for  high  values  of  0 
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Figure  3.4:  The  temperature  distribution  across  a  class  of 

stationary  slabs  with  uniform  composition  is 
given.  Here  the  stress  strain-rate  relations  is 
given  by  equation  (3.7;  and  the  temperature  is 
given  in  term  of  the  dimensionless  temperature 
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3. 4  Application  to  the  Seismic  Low  Velocity  Zone 

The  velocity  of  seismic  waves  in  general  increases  with 
pressure  and  decreases  with  increasing  temperature.  One  or 
the  other  of  these  two  effects  can  dominate  for  certain 
temperature  and  pressure  gradients  depending  on  the 
material.  In  fact  many  early  explanations  of  the  existence 
of  the  seismic  low  velocity  zone  were  soley  dependent  on 
this  phenomena  (Birch  p.  1952,  Valle  P.E.  1956,  MacDonald 
G.J.F,  and  Ness  N.F.  1961)  .  The  critical  temperature 
gradients,  particularly  for  S  waves  are  probably  exceeded 
within  the  outer  150  km  so  there  is  little  doubt  that  this 
explanation  is  at  least  partially  correct  (Ringwood  LE« 
1975,  Lieberman  R.C.  and  Schrieber  E.  1969).  It  is  however 
pointed  out  by  Ringwood  (1975)  that  chemical  and 
mineralogical  heterogeneity  must  also  play  a  part  since  the 
largest  vertical  temperature  gradients  in  the  earth  are 
immediately  below  the  Mohorovic  discontinuity  where  the 
seismic  velocity  still  increases  with  depth.  I  will  not 
consider  these  factors  here  as  they  require  a  rather 
detailed  petrological  study  which  is  given  by  Ringwood 
(1975)  . 

The  general  physical  situation  which  results  in  plate 
formation  is  extremely  complex  and  open  to  much  dispute.  The 
theories  in  this  area  vary  from  the  crust  and  lithosphere 
being  formed  by  chemical  separation  (Ringwood  A.E.  1975)  to 
a  fairly  uniform  composition  throughout  with  the  vertical 
inhoraogenieties  being  caused  by  phase  changes  resulting  from 
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the  increasing  pressure  with  depth  (Kennedy  G.C.  1959) .  Here 
I  simply  mention  the  situation  as  I  envisage  it  from  a  fluid 
dynamics  point  of  view.  In  the  formation  of  the  plates  there 
are  a  number  of  factors  which  must  be  considered  such  as  the 
melting  point,  conductivity  and  density  of  the  different 
materials  involved.  In  general  the  least  dense  material,  the 
material  of  lowest  melting  point  and  the  material  of  lowest 
conductivity  all  tend  to  rise  to  the  top  in  lithosphere 
differentiation.  The  crust  is  then  formed  by  material  with 
some  or  all  of  these  properties  depending  on  the  degree  to 
which  any  of  these  properties  can  dominate.  The  rest  of  the 
lithosphere  is  then  formed  mainly  by  material  with  the 
opposite  properties  which  has  been  depleted  of  the  crustal 
material.  Once  the  plates  have  been  formed  they  separate  the 
high  temperatures  of  the  earth's  interior  from  the  much 
cooler  temperatures  at  the  surface.  This  layer  of  thermal 
insulation  is  very  efficient  because  its  high  viscosity 
allows  heat  to  escape  only  by  conduction  over  the  majority 
of  the  earth's  surface. 

Below  the  plates  convection  of  the  earth's  mantle  then 
causes  smaller  temperature  gradients.  This  description 
however  is  complicated  by  the  fact  that  the  thickness  of  the 
plates  varies  from  less  than  50  km  near  the  ridges  to  250  km 
within  some  of  the  continents  (Isacks  B.  et  al  1968) .  This 
causes  the  temperatures  in  the  continental  lithosphere  which 
are  at  comparable  depths  to  the  seismic  low  velocity  zone 
beneath  the  oceans  to  be  approximately  200°  cooler.  As  a 
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result  the  marginal  stability  beneath  the  oceans  would  not 
be  experienced  in  the  mantle  material  beneath  the  continents 
in  an  unstressed  state. 

The  pressure  in  the  seismic  low  velocity  zone  beneath 
the  continental  lithosphere  is  in  fact  very  similar  to  that 
of  the  mantle  material  beneath  the  seismic  low  velocity  zone 
under  the  oceans.  One  thus  does  not  get  the  same  condition 
of  an  ambient  temperature  close  to  the  melting  temperature 
because  of  the  effect  of  pressure  on  increasing  the  melting 
temperature.  A  region  of  marginal  stability  therefore  does 
not  exist  to  the  same  degree  as  beneath  the  oceans  allowing 
a  sharp  region  of  decoupling  of  motion  of  the  plates  and  the 
convection  cells.  As  a  result  the  seismic  low  velocity  zone 
beneath  the  continents  is  far  more  complex  than  beneath  the 
oceans, 

\ 

I  do  not  propose  to  study  this  region  beneath  the 
continents  in  depthf  but  wish  to  simply  consider  some  of  its 
properties.  From  the  present  knowledge  of  plate  motions 
relative  to  the  mantle  convection  cells  it  is  generally 
acknowledged  that  the  plates  which  contain  continents  move 
slower  than  those  without.  This  may  be  concluded  by  looking 
at  relative  plate  motions  and  considering  the  hot  spots 
which  are  assosiated  wh it h  uplift  (and  thus  can  be 
considered  as  originating  from  depth)  as  being  essentialy 
fixed  ^Burke  K.  et.  al.  1973).  This  observation  then  seems 
to  imply  that  the  continents  impart  greater  amount  of  shear 
stress  to  the  ast henosphere  and  thus  impede  plate  motions 
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But  partial  melting  is  assumed  to  be  the  dominant  cause 
of  the  low  viscosities  in  the  seismic  low  velocity  zone 
beneath  the  oceans  (this  has  also  been  speculated  by  many 
others  e.g.  Cathles  1. M.  1975).  Thus  in  order  to  conclude 
that  the  shear  stress  in  the  seismic  low  velocity  zone 
beneath  the  continents  is  greater  one  must  also  concJ  ude 
that  the  degree  of  partial  melting  is  also  less.  This 
follows  directly  from  the  fact  that  one  must  have  larger 
viscosities.  This  however  is  about  as  far  as  I  can  argue  for 
the  continents  in  general.  Due  to  the  tremendous  differences 
of  the  tectonic  environments  of  the  different  continents  one 
finds  that  for  any  general  theory  developed  for  the 
existence  of  the  seismic  low  velocity  zone,  almost  every 
individual  continent  is  an  exception.  This  is  a  result  of 
such  phenomena  as  glacial  rebound,  continental  collisions, 
the  subduction  of  oceanic  plates  beneath  the  continents, 
large  variations  in  the  depth  of  different  continents  etc. 
all  of  which  have  substantial  changes  on  the  physical 
conditions  in  the  seismic  low  velocity  zone. 

Consider  now  the  lithosphere  -  asthenosphere  boundary 
beneath  the  ocean.  Here  one  has  a  slight  decrease  in  density 
from  the  lithosphere  above  (of  about  3.35  gm/cm3  to  3.3 
gm/cm3)  and  much  larger  temperature  gradients  than  beneath 
the  continents.  This  causes  the  material  to  be  subjected  to 
relatively  high  temperatures  which  are  close  to  the  melting 
point  of  the  material  at  the  pressures  at  those  depths.  This 
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material  may  then  be  considered  as  marginally  stable  meaning 
that  it  is  either  close  to  partial  melting  or  slightly 
beyond  the  solidus  of  some  of  the  components. 

The  results  of  this  section  will  indicate  that  the 
thermo-mechanical  coupling  of  the  plate  motions  above  to  the 
more  viscous  mantle  below  cannot  be  solely  responsible  for 
the  observed  low  viscosities  in  the  seismic  low  velocity 
zone.  However  when  an  increasing  viscosity  with  depth  due  to 
changing  physical  properties  is  taken  into  consideration  the 
shear  heating  can  then  result  in  the  observed  changes  in 
size  and  intensity  of  the  seismic  low  velocity  zone.  This 
can  occur  due  to  changes  in  the  relative  motion  of  the 
plates  and  mantle  convection  cells.  Also  it  can  be  seen  that 
shear  heating  must  at  least  be  responsible  for  an  increase 
in  the  amount  of  melt  in  existence  in  the  seismic  low 
velocity  zone.  This  melting  may  also  be  aided  by  the 
presence  of  up  to  .1%  water  (hingwood  A.  E.  1975}.  The 
increasing  viscosity  with  depth  is  then  a  result  of  the 
pressure  gradients  becoming  greater  than  the  temperature 
gradients  and  a  changing  composition  due  to  a  rise  of  the 
melting  material.  I  will  show  at  the  end  of  this  section 
that  one  may  construct  stationary  solutions  for  the  material 
becoming  more  viscous  with  depth  due  to  compositional  and 
pressure  changes.  The  molten  material  generated  in  the 
seismic  low  velocity  zone  then  rises  and  cools  as  part  of 
the  lithosphere.  This  fits  theories  of  lithosphere 
thickening  (e.g.  Parker  R.  L.  and  Oldenburg  D.  W.  1  973)  . 
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The  seismic  low  velocity  zone  as  a  region  of  steady 
state  thermo- mechanical  coupling  of  a  marginally  stable 
material  would  at  first  sight  seem  to  imply  that  the  low 
velocity  zone  may  be  driven  to  instability  through  the 
application  of  a  small  amount  of  shear  stress  by  the  moving 
plates.  This  however  is  certainly  not  the  case.  Consider  the 
Bridgman  (1935,  1936,  193?)  experiments  once  again.  They 
illustrate  that  materials  with  low  melting  points  deform 
smoothly  without  experiencing  a  degree  of  thermo-mechanical 
coupling  which  leads  to  dramatic  stress  relief.  A  thermal 
instability  is  thus  just  localized  runaway  heating  which 
occurs  in  an  ever  narrowing  region  due  to  the  shear  heating 
increasing  faster  than  the  ability  of  the  material  to 
conduct  heat.  Upon  examination  however  one  finds  that  the 
low  velocity  zone  must  be  even  more  stable  than  a  material 
with  low  melting  point.  This  is  because  it  is  not  composed 
of  one  material  close  to  melting  but  numerous  materials 
(olivine,  pyroxine,  garnet,  and  perhaps  in  some  regions 
amphiboie  (Ringwood  A.E.  1975))  which  melt  over  a  range  of 
temperature.  Thus  as  the  solidus  for  some  of  the  materials 
is  passed  the  viscosity  drops  resulting  in  a  smaller  amount 
of  shear  stress  being  imparted  by  the  motion  of  the  plates. 


This  is  simply  a  consequence  of  eqn  (3.3)  or  (3,9)  The 
majority  of  the  material  then  remains  unmelted  with  the  melt 
being  considered  at  present  as  uniformly  distributed  across 
the  low  velocity  zone.  The  stress  being  imparted  thus  is 
able  to  be  transported  at  large  velocities  across  the  slab 
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as  a  vrhoie  making  the  concept  of  a  localized  instability 
even  more  unrealistic* 

One  possible  generalization  of  the  work  in  this  chapter 
now  becomes  obvious.  That  would  be  to  introduce  the  concept 
of  a  concentration  field  of  melt  which  would  allow  the 
material  to  reach  the  solidus  over  a  range  of  temperatures. 
The  equations  for  this  case  are  shown  in  appendix  B.  The 
problem  in  this  case  is  that  one  must  consider  the  behaviour 
of  three  coupled  equations  which  makes  finding  solutions 
very  difficult. 

In  order  to  evaluate  the  results  of  this  chapter  the 
following  values  which  appear  to  be  representative  of 
seismic  low  velocity  zone  conditons  are  used,  f Temperatures 
are  in  °Ce  other  units  are  cgs  except  where  noted.) 

c  =  U  x  10?  y0  -  1022 

\. 

TB  =  1100  E  =  60kcal/mole 

Tm  =  1400  E/E  =  3  x  10* 

P  =  3.5  R  =  2  x  10  kcal/°K/ifiole 

k  =  3  x  10s  n  -  3 

I  now  calculate  the  shear  stress  necessary  to  maintain  this 
shear  coupling  zone  in  a  marginally  stable  state  over  a 
distance  of  100  km»  This  zone  is  assumed  to  be  caused  by  the 
shear  stress  applied  by  the  moving  plate  at  the  top  and  the 
more  viscous  material  below. 

Using  the  stress  strain  rate  relation  (2.2)  the  stress 
in  this  simple  stationary  thermo-mechanical  system  is  given 
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by 
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(3.13) 


If  the  temperature  is  a  monotonical ly  increasing  function 
with  depth  then  £  is  the  thickness  of  the  low  velocity 
zone  and  yQ  is  its  viscosity  at  the  plate  boundary.  Egn. 
(3.13)  then  yields  an  applied  stress  of  approximately  50 
bars  which  in  turn  yields  a  strain  rate  at  any  position 
across  the  slab  of 


e(x)  =  5  x  10  1j  exp ( $ (x) ) /sec 


(3.14) 


Here  p  (x)  goes  from  0  to  1.2  which  corresponds  to  a 
temperature  range  of  1100°C  to  1220°C. 

If  one  now  assumes  a  plate  motion  of  5  cm/year,  then  a 
velocity  of  approximately  IQ-7  cm/sec  is  obtained  for  the 
top  of  the  low  velocity  zone.  Also  if  a  much  smaller 
velocity  is  assumed  for  the  more  viscous  material  below  and 
a  constant  velocity  gradient  across  the  low  velocity  zone 
then  a  strain  rate  of  approximately  10-**  /sec  results.  This 
result*  is  almost  identical  to  the  strain-rate  obtained  in 
eqn(3.l4)  and  depending  only  on  the  physical  constants  and 
the  assumption  of  marginal  stability  for  the  equations  in 
section  (3.2)  . 

Now  in  order  to  use  the  stress  strain  rate  relation 
(3.7)  one  must  first  know  and  f>  c  so  that  the  stress  may 
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be  calculated  from  eqn  (3.10).  Here  a  temperature  range  of 
1100°C  to  1400°C  is  taken  which  is  probably  more  reasonable 
than  the  range  obtained  in  the  previous  example.  In  this 
case  the  stress  once  again  works  out  to  be  approximately  50 
bars  and  the  strain  rate  is  approximately 
5  x  10“ 15  x  exp  (9. 1-1/0  )  /sec 

Here  1/$>  goes  from  7.1  to  9.1  and  thus  the  strain-rate  once 
again  fits  very  well  with  observed  plate  motions. 

An  increasing  pressure  and  the  rise  of  molten  material 
both  essentially  cause  the  value  of  y  Q  to  be  an  increasing 
function  of  depth.  Thus  I  wish  to  consider  the  behavior  of 
such  a  situation  in  terms  of  a  value  yQ  (x) ,  for  which  the 
previous  equation  may  still  be  solved  analytically. 

Assume 

V  =  yQ  exp  (-(j))  exp  [{1-0x1 


where  £  =  0  is  the  bottom  of  the  low  velocity  zone,  y  =102i 

and  x  is  a  constant  such  that  yG  x  exp(X)=102?.  We  may 
then  write  eqn  (3.4)  in  the  form 
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-A  exp  (<{>-  (1-0  x) 


How  let  y  ~<p  -(1-  K  )X  and  we  then  have 


which  has  the  solution 
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$  =  -xC-2  In  /p- 
Here 


cosh  <  cosh 


-1  Z 
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and  3  $  /9£  |  is  the  temperature  gradient  evaluated  at  the 
boundary  of  the  slab.  Thus 


M. 


b 


=  -X  +  E 2  tanh  < 


cosh 


exp 


X 

2 


From  the  value  of  p  at  £  =1  we  can  solve  for 
a<f>  /9£  |  .  We  may  then  determine  the  position  of  the 
maximum  temperature  by  setting  8  p  /9£  ~0  and  solving  f or  £  . 
The  effect  of  this  assumption  of  a  composition  gradient  is 
that  it  causes  the  maximum  temperature  to  be  shifted  towards 
the  boundary  of  lower  viscosity  (see  figure  (2.5)).  With  the 
introduction  of  heating  from  below  (the  boundary. of  higher 
viscosity)  one  may  get  rid  of  this  maximum  temperature  as 
illustrated  in  figure  (2.  7)  . 
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Chapter  4 


The  Inertial  Field  Equations  and  Their  Application  to 
Deep  Earthquakes 

4xi  A  Physical  Discussion 

Studies  of  intermediate  and  deep  earthquake  mechanisms 
indicate  that  the  principal  stresses  are  developed  within 
the  down  dipping  plate  and  in  the  direction  of  plate  motion. 
(Sykes  L 0  R *  1966,  Isacks  B.L.  and  Molnar  P.  1969,  1971  , 
Isacks  B.L.  et  al  1969)  .  However  the  proposed  physical 
explanations  of  deep  earthquake  mechanisms  all  have 
problems.  Consider  first  a  proposal  of  Isacks  and  Koinar 
(1971) £ 

If  one  assumes  that  the  slab  cannot  penetrate  the  650 
km  discontinuity  then  the  upward  transmission  of  the 
compressional  stress  could-  cause  deep  earthquakes.  This 
mechanism  is  rejected  in  great  detail  by  Ringwood  (1975). 
Also  it  has  been  shown  by  Griggs  (1972)  that  such  resistance 
would  result  in  buckling  of  the  subducting  plate,  for  which 
there  is  little  evidence,  and  thus  could  not  transmit  the 
necessary  stress  for  hundreds  of  kilometers  upwards  in  any 
case. 

Another  proposal  has  been  made  by  Ringwood  (1975)  and 
others  (Bridgman  P.W.  1945,  Ringwood  A.E.  1956,  1967,  Evison 
F.  1963,  1967,  Dennis  I.C.  and  Walker  C.  1965).  Here 
intermediate  and  deep  focus  earthquakes  are  assumed  to  be 
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caused  by  phase  transformations  associated  with  sudden 
volume  changes  in  the  subducting  slab,  Two  problems  with 
this  mechanism  have  been  suggested  by  Ringwood  (1975).  It  is 
not  understood  how  phase  transitions  can  occur  quickly 
enough  to  cause  an  earthquake  and  the  analyses  of  first 
motions  from  earthquake  seismograms  show  douple  couple 
mechanisms  exactly  as  in  shear  failure  rather  than  monopolar 
implosion  as  has  been  believed  to  be  required  in  phase 
change  mechanisms.  Ringwood  (1975)  then  simply  passes  the 
first  problem  off  as  due  to  the  large  temperature  gradients 
and  the  second  problem  as  the  relief  of  pre-existing  stress 
due  to  gravitational  body  forces.  An  analysis  of  the  seismic 
moments  of  intermediate  and  deep-forcus  earthquakes  has  been 
done  by  McGarr  (1977).  He  claims  to  get  reasonable 
consistency  under  the  assumption  of  phase  change  induced 
earthquakes  in  four  island  arc  regions. 


I  am  however  not  aware  of  any  research  which  has  been 
done  on  the  dynamics  of  this  mechanism.  Part  of  the  reason 
for  this  apparent  absence  of  work  in  this  area  is  certainly 
due  to  the  fact  that  the  theoretical  base  consists  in 
applying  a  static  theory  to  a  dynamic  system.  I  would  like 
to  suggest  that  with  our  present  meagre  knowledge  of  the 
fashion  in  which  the  phase  changes  occur  it  is  just  as 
reasonable  to  conclude  that  they  inhibit  deep  earthquakes  as 
cause  them.  I  make  this  suggestion  since  a  viscosity 
increase  is  probably  associated  with  the  majority  of  the 
phase  changes  in  the  subducting  slab  because  of  the  density 
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increases  which  occur  and  I  am  aware  of  no  evidence  that  the 
phase  changes  can  occur  as  quickly  as  necessary  to  induce  an 
earthquake  under  the  conditions  present  in  the  subducting 
slab.  Also,  almost  every  subducting  plate  on  earth  exhibits 
an  absence  of  earthquakes  in  the  region  where  the  Olivine- 
Spinel  phase  change  is  thought  to  occur  in  the  slab. 

It  has  also  been  suggested  that  deep  earthquakes  may  be 
evidence  of  a  fluid  instability  within  the  subducting  slab, 
caused  by  shear  stress  imparted  by  the  surrounding  mantle 
(Griggs  D.T.  and  Baker  D.U.  1969,  Shaw  H.K,  1989).  Although 
this  theory  would  seem  to  fit  best  with  observations  of 
focal  mechanism  of  deep  earthquakes  it  is  shown  in  sec. (4.3) 
that  the  subducting  slab  is  orders  of  magnitude  away  from 
instability. 

I  wish  to  look  further  into  the  possibility  of  fluid 
instability  as  a  mechanism  for  deep  earthquakes  in  this 
chapter.  First  one  should  consider  carefully  where  the 
instability  would  occur.  I  show  in  section  (4.3)  that  the 
only  region  which  may  be  realistically  driven  to  instability 
is  the  boundary  layer  between  the  mantle  and  upper  edge  the 
subducting  slab.  The  interior  of  the  slab  is  far  too  stable 
to  exhibit  instability  while  the  surrounding  mantle  is  able 
to  deform  smoothly.  The  boundary  layer  on  the  other  hand  has 
the  greatest  amount  of  heat  production  (Turcotte  D.  L.  and 
Schubert  G.  1973  also  see  sec.  4.4)  and  is  initially  far 
from  melting.  These  are  the  precise  conditions  necessary  for 
the  generation  of  a  thermal  instability  according  to  the 


to  w*  *  •  ' 


s 


J  .  4  *  ■ 


80 


experimental,  observations  made  by  Bridgman, 

4, 2  An  Analysis  of  the  Equations  of  Motion 

I  consider  here  a  discription  of  the  inertial  effects 
of  a  highly  viscous  fluid  when  subjected  to  a  time  dependent 
boundary  stress.  I  start  by  discussing  a  fluid  contained 
within  an  infinite  horizontal  slab  of  depth  21  »  The 
analysis  I  use  for  the  propagation  of  the  velocity  field  in 
a  highly  viscous  fluid  is  similar  to  the  theory  for  the 
propagation  of  a  plastic  wave  in  a  visco-elastic  solid 
(Taylor  G.  I.  1946,  von  Karman  ?.  and  Duwez  P.  1950, 
Rakhmatuiin  K.  A.  1945).  In  the  case  of  a  visco-elastic 
solid  one  is  dealing  with  a  solid  with  some  fluid 
properties,  i.  e.  the  "fluid  like”  propagation  of  a 
displacement.  The  difference  here  is  that  I  am  dealing  with 
a  fluid  with  some  solid  properties,  i.e.  the  "solid  like" 
propagation  of  a  velocity.  Both  analyses  are  identical  in 
that  they  represent  the  propagation  of  a  stress  wave  across 
the  material.  The  difference  is  simply  the  method  by  which 
stress  is  propagated  in  the  two  cases. 

If  a  uniform  shear  stress  is  applied  at  the  boundaries, 
that  stress  will  be  assumed  to  propagate  across  the  slab  at 
the  velocity  3(x,t).  Here  x  is  the  vertical  co-ordinate  and 
t  is  the  time.  8(x,t)  therefore  represents  the  ratio  A x/f  t 
of  the  distance  a  velocity  perturbation  will  propagate, 
orthogonal  to  the  direction  of  the  perturbation,  to  the  time 
At  of  propagation  (i. e.  in  the  case  of  a  fluid  of  constant 
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viscosity  throughout  and  initially  constant  velocity 
throughout  v  (x<A  x ,  t+A  t )  -  v  {x,  t)  )  .  As  each  part  of  the  slab 
then  becomes  subjected  to  a  horizontal  force  we  obtain  a 
corresponding  motion.  The  slab  thus  undergoes  plastic 
deformation  and  we  obtain  a  continuous  distribution  of  heat 
sources.  This  continuum  of  heat  sources  is  the  shear  heating 
due  to  the  viscous  dissipation  of  the  forces  responsible  for 
the  dynamics.  The  heating  of  the  fluid  then  causes  a 
decrease  in  the  viscosity  and  thus  the  ability  of  the  fluid 
to  dissipate  the  shear  stress  imparted  to  it  (i.e.  the 
plastic  shear  wave  velocity  decreases)  .  This  causes  a  change 
in  the  dynamics  which  according  to  the  previous  argument 
once  again  changes  the  thermal  field.  One  thus  obtains  a 
coupled  thermo- mechanical  system  in  which  both  fields  must 
be  solved  simult aneously  except  for  the  special  case  of 
stationary  solutions,  i.e.  solutions  which  do  not  contain 
any  inertial  motions,  which  were  discussed  in  Chapter  3. 

Now  consider  the  evolution  of  an  instability  in  such  a 
material.  In  Appendix  C  it  is  shown  that  one  cannot  obtain  a 
stable  stationary  solution  for  a  homogenous  half  space 
subjected  to  constant  shear  no  matter  how  slowly  the  shear 
stress  is  increased  to  the  desired  amount.  One  should 
however  remember  that  this  is  strictly  a  mathematical 
result.  When  using  the  idealized  case  of  a  homogeneous  half 
space,  the  above  result  can  only  be  meaningful  where  the 
distance  from  the  boundary  to  the  position  of  the 
instability  is  a  distance  over  which  the  basic  physical 
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assumption  can  be  considered  reasonable  and  the  time  for  an 
instability  to  evolve  is  less  than  the  time  for  the 
phenomena  being  modeled  to  occur. 

It  should  be  noted  that  all  models  of  deep  earthquake 
mechanisms  are  very  speculative  due  to  the  great 
uncertainties  in  knowledge  of  the  geologic  and  physical 
conditions.  Here  I  attempt  to  find  a  simple  model  which  fits 
with  the  present  experimental  evidence  and  reguires  a 
minimum  number  of  hypotheses.  I  wish  to  consider  the 
equations  of  fluid  flow  in  more  detail  in  the  remainder  of 
this  section  and  in  section  (,4.3)  before  presenting  a  new 
model  for  deep  earthquakes  in  section  (4.4)  based  on  these 
equations. 

The  equation  of  motion  of  a  horizontally  sheared 
viscous  fluid  in  the  absence  of  convection  is  given  by 


9  (pv)  _  9tt 

9 1  ~  9x 


(4.1) 


Here  the  velocity ,  vf  is  in  the  y  direction  horizontal  to 
the  surfaces,  x  is  perpendicular  to  the  surface,  p  is  the 
density  and  II  is  the  momentum  flux  density.  In  the  case  of 
pure  shear  (i.e.  when  the  momentum  flux  due  to  convection 
and  pressure  gradients  is  zero)  one  may  write 

it  =  -a 


where  o  is  the  shear  stress  and  -a  represents  the  momentum 
flux  density  due  to  viscous  dissipation  for  pure  shear. 
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In  the  absence  of  phase  diffusion  the  equation  of  heat 
transfer  is  given  by 


pc 


9T 
9 1 


a  +  |_ 
dx  9x 


jr  9T 

K  3^ 

L. 


(4.2) 


The  right  hand  side  of  this  equation  may  also  be  written  in 
the  form 


Applying  Gauss’s  theorem  to  a  unit  volume  of  the  fluid  we 
find  k  0T/8x  represents  the  energy  flux  density  through  the 
surface  due  to  thermal  conduction;  va  represents  the  energy 
flux  density  through  the  surface  due  to  non- inertial  viscous 
effects;  and  v9a  /9x  represents  the  net  energy  flux  density 
in  the  unit  volume  due  to  the  inertial  coupling  of  the 
thermal  and  mechanical  field  equations.  The  last  term  may 
now  be  divided  as  follows 


v 


9a 

9x 


pv 


9  v 

9 1 


=  P 


3 (x, t) v 


9  v 
9x 


( P  3 (x, t) v2) 


p  v2 
9x 


(4.2a) 


where  the  time  of  propagation  for  the  dynamics  resulting 
from  an  applied  force  is 
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3 (xf t) 


dx 


(4.3) 


and  3(x,t)  is  the  velocity  of  propogation  of  the  dynamics 
in  the  direction  orthogonal  to  the  applied  force.  Here 
p3ix„t)v2  represents  the  energy  flux  density  within  a  unit 
volume  of  the  fluid  due  to  the  inertial  interactions  of  the 
thermal  and  mechanical  field  equations  in  that  volume. 

Assuming  gT/St-O  the  net  energy  flux  in  a  unit  volume 
of  the  fluid  becomes  zero.  This  may  be  written  as 


6 {X/t)  37  =  0  (4.4) 

since  the  inertial  heating  is  strictly  a  result  of  a  heat 
source  due  to  the  dynamics  and  thus  either  $(x,t)=0  or 
3T/3x-0.  If  $(xrt)=0  then  one  can  see  from  (4.2a) 
immediately  that  8v/8t=0.  If  8T/3x=0  then  the  net  heat 
conduction  must  be  zero  and  the  heat  production  must  be  the 
same  everywhere.  Now  assuming  a  Newtonian  stress  strain-rate 
relation  we  have  y  { 9v/ 3x) 2=constant  and  since  y  is  a 
function  of  ?  only  y-constant.  Thus  3  v/3x~  constant  and 
since  the  heat  production  may  also  be  written  as 

o3v/3x~constan t,  we  have  cr=constant.  Thus 

9v  _  1  3o 
3t  p  3x  “  0 


Now  assume 


3v/3t=0«  Differentiating  o  =y  3v/3x  with  respect 
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to  t  we  obtain 

9cT  _  _9jJ  9v 
9 t  9 1  9x 

which  may  be  written  in  the  form 
OR  (X  t)  lX  =  IX  1±  IX 

'  '  gt  ^  gt  8x 

Thus  we  have 

n  =  lX  li5  lX 
U  “  3<J)  9 1  9x 

but  9^/9$  ^0  since  viscosity  is  a  strong  function  of 
temperature  in  a  highly  viscous  fluid  and  9v/9x^O  since  we 
have  a  viscous  fluid.  Thus  any  applied  stress  must  result  in 
velocity  gradients  across  the  slab.  Therefore  for  a  highly 
viscous  fluid  undergoing  shear  stress  one  obtains  the 
condition 


li. 

9 1 


0 


lx 

9 1 


0 


We  thus  have  the  result  that  heating  can  occur  only  if  there 
are  accelerations  and  accelerations  can  occur  only  if  there 
is  heating. 


4^3  The  Transition  to  the  Inertial  Field  Equations. 


Consider  now  the  solutions  which  evolve  for 

horizontal  slab  of  half -t hickness  l  when  a 

n  (t)  is  applied  at  each  boundary.  The  forces 
b 

for  the  stresses  are  assumed  to  be  in  opposite 


an  infinite 
shear  stress 
responsible 
directions. 
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The  rate  of  heat  production  due  to  viscous  shear  heating  at 
any  point  across  the  slab  is  once  again  given  by  oi  where 
a  and  £  are  the  stress  and  strain  rate  at  that  point 
respectively. 


A  stress  strain-rate  relation  which  holds  for  most 
geologic  material  in  the  neighbourhood  of  instability  was 
given  by  equation  (3.7).  However  I  first  wish  to  consider 
the  simpler  relationship  obtained  for  T  =Tb  in  eqn  {2,2}  and 
the  newtonian  approximation  eqn  (2.2a). 


The  heat  balance  equation  is  then  given  by 


9<J) 

at1 


o 


(£  ,  t ' )  exp  (<j) ) 


where 

f  =  a(T-T0')  '  C  =  x/£  and  t'  = 

'  '  pcH2 

Momentum  balance  is  given  by  the  equation 

9v  _  pc~  9 a  ( g  f t  * ) 

9 1 ' 

Here  v  is  the  dimensionless  velocity. 

In  the  absence  of  inertial  term  it  was  shown  in  chapter 
3  that  analytic  solutions  may  be  obtained  for  A<.88.  Now  in 
order  to  illustrate  the  importance  of  the  inertial  terms,, 
the  evolution  of  the  system  from  its  initial  condition  to 
melting  is  followed  by  different  processes.  First  consider  a 
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quasi -stat ionary  process*  This  implies  that  we  study  the 
evolution  of  the  system  for  an  infinitely  slow  rate  of 
increase  of  stress  at  the  boundary.  The  solutions  at  any 
time  are  then  given  by  the  stationary  theory.  Assume  the 
equations  yield  solutions  for  values  of  $> r  ^ M  tfhere  $c 
is  the  central  temperature  and  $>M  is  the  melting 
temperature.  The  maximum  temperature  occurs  at  the  centre  of 
the  slab  and  it  is  there  that  the  condition  $=<PM  is  first 
encountered  and  the  slab  begins  to  melt.  This  onset  of 
melting  and  its  attendant  reduction  in  viscosity  creates  an 
instability  for  AC. 88. 

Now  consider  the  same  evolution  as  before  with  0  c  <  <?>  M 
for  all  possible  stationary  solutions.  The  position  of  the 
instability  may  be  found  by  considering  the  slab  as  three 
separate  slabs 

(Slab  1:  -1^  5^-e  ,  Slab  2:  -e  ^  f  Slab  3:  e  S^l) 

Here  e  is  an  arbitrarily  small  number.  If  A  =.88  for  the 
slab  as  a  whole,  slabs  1  and  3  once  again  yield  stationary 
solutions  (see  appendix  C)  capable  of  transporting  a  greater 
amount  of  stress.  Slab  2  has  a  central  maximum  temperature 
of  <§>  =1.2  which  must  correspond  to  A=.88  in  this  slab  (see 
sec.  3.3).  Thus  any  stress  increase  in  this  slab  must  result 
in  an  instability.  By  this  analysis  one  has  now  obtained  the 
result  that  any  slab  which  is  composed  of  one  material  and 
whose  motion  and  heating  is  determined  by  equation  (2.2)  can 
reach  instability  only  at  the  center.  The  generalization  of 
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this  argument  to  the  more  complex  stress  strain  -  rate 
relation  (3.7)  is  trivial  since  stationary  solutions  exist 
for  any  value  of  the  applied  stress  (sec  3,4)  and 
instability  simply  occurs  when  the  central  maximum 
temperature  reaches  the  melting  temperature. 

Note  that  here  one  is  assuming  that  one  has  a  perfect3.y 
homogeneous  fluid  composed  of  a  single  material  and  a  very 
well  defined  melting  point.  Here  one  is  also  considering  a 
constant  stress  boundary  condition  and  not  a  constant 
velocity  boundary  condition  in  which  case  the  stress  can 
usually  evolve  to  a  sufficiently  low  value  that  instability 
will  not  occur.  This  analysis  thus  does  not  in  any  way 
contradict  the  work  in  section  (3.4). 

Now  consider  the  instantaneous  application  of  stress  at 
the  boundary.  Here  a  b  has  a  time  dependence  like  a  H (t) 
where 

H  (t)  =[  1  t>l 
0  t <0 

and  a  is  a  constant.  Since  the  stress  must  propagate  into 
the  slab  at  a  finite  velocity  8  se  have  the  condition 
3a  /  3x  |  =  5 ( |  x !  )  where  S(x)  is  the  dirac  delta 

function.  This  implies  that  v  has  a  space  dependence  like 
a  h ( | x (  -  £  )  and  3  v/ 3  x  has  a  space  dependence  like 
6  ~Z  )  for  -£  <x<  Z  at  time  zero.  Since  the  heat 

production  term  is  the  product  of  a  finite  stress  and  an 
infinite  strain  rate,  the  density  of  heat  sources  at  the 
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boundary  is  infinite. 

With  an  infinite  density  of  heat  sources  the  viscosity 
drops  infinitely  rapidly.  This  means  the  stress  cannot 
propagate  into  the  slab  and  the  velocity  gradient  must 
increase  without  limit.  In  particular,  melting  will  now 
occur  on  the  boundary  and  not  at  the  centre.  A  true  physical 
system  is  neither  quasi-stationa  ry  nor  has  a  zero  time 
interval  for  the  application  of  stress.  With  a  finite  speed 
of  applciation  of  stress  instabilities  can  be  expected  at 
any  distance  between  the  boundaries  and  the  center  depending 
on  the  rate  of  increase  of  stress  and  the  magnitude  of  the 
stress.  This  indicates  that  the  characteristic  distance  to 
instability  is  a  function  cf  the  time  rate  of  change  of  the 
boundary  force,  the  material  properties  and  the  initial 
conditions. 

v 

I  wish  to  consider  now  the  physical  situations  which 
lie  between  these  two  very  unphysical  extremes.  Consider  an 
infinite  slab  of  thickness  Z  which  is  held  at  a  constant 
temperature  at  one  surface  and  thermally  insulated  at  the 
other  surface.  The  thermally  insulated  surface  thus 
corresponds  to  the  centre  of  the  slabs  discussed  previously. 
Assume  the  slab  is  subjected  to  an  applied  force  f  (t) 
horizontal  to  the  cooled  surface.  The  traction  at  the 
boundary  is  given  by 

ab  (t)  *n  =  ?  (t) 

Here  the  only  non  zero  component  of  traction  yields  a 
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horizontal  shear  stress.  This  stress  propagates  to  any  depth 
x  in  the  slab  in  a  time 


t 

P 


(V 


X 


o 


3  (T  (x  ,  t )  ) 


dx 


Thus  we  may  write 

cr(x,t)  =  a^(t-t  (x)) 


where  a(xrt)  represents  the  shear  stress  at  the  position  x 
within  the  slab  at  time  t.  At  this  point  it  should  be  noted 
that  one  may  obtain  further  inertial  effects,  such  as  those 
which  can  result  from  the  partial  attenuation  of  the  stress 
waves,  which  have  not  been  considered*  Neglecting  these 
effects  should  however  give  a  reasonable  upper  bound  to  the 
depth  and  time  to  instability. 


The  equations  of  motion  may  now  be  written  in  terms  of 
the  boundary  stress  and  shear  wave  velocity  as 
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(4.5) 
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The  velocity  of  propagation  of  the  dynamics,  $(xft) 
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will  here  be  considered  to  be  the  shear  wave  velocity  when 
8a  /  3 t  is  large  since  in  this  case  the  material  should 
exhibit  only  solid  properties  and  the  plastic  wave  velocity 
when  do/  9 t  is  small  since  the  material  should  exhibit 
only  fluid  properties.  It  is  however  such  a  highly  nonlinear 
function  of  temperature  (especially  in  the  neighbourhood  of 
melting)  that  it  would  be  foolish  to  try  to  solve  eg*s  (4.5 
and  4.6)  in  any  generality  either  analytically  or 
numer  ically. 

One  case  however  which  may  be  calculated,  is  the  time 
for  the  stress  to  propagate  to  any  point  in  the  interior 
when  a  stationary  solution  is  subjected  to  a  small 
perturbation  in  the  applied  stress  at  the  boundary.  This 
case  represents  a  reasonable  starting  point  for  considering 
the  onset  of  instability.  Since  the  viscosity  decreases 
exponentially  with  increasing  temperature  it  is  reasonable 
to  assume  the  velocity  of  propagation  of  the  dynamics  also 
decreases  exponentially  with  increasing  temperature. 
Viscosity  is  after  all  simply  a  measurement  of  the  ability 
of  a  fluid  to  diffuse  momentum.  Therefore  with  the 
approximation  T  -  TB  and  the  Newtonian  stress  strain-rate 
relation  one  obtains  a  time 


for  a  small  perturbation  in  the  boundary  stress  to  propagate 
to  a  distance  £  from  the  thermally  insulated  boundary. 
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Kere  the  shear  wave  velocity  is  assumed  to  be  of  the  form 

3  ( 4> )  =  3q  exp  (-(f)) 

where  is  the  velocity  of  propagation  of  the  dynamics  at 

<£>  =  0.  This  form  for  $  does  not  allow  the  stress  to 
propagate  past  an  instability.  It  also  fits  with  the  natural 
expectation  that  3  have  the  same  temperature  dependence  as 
viscosity.  The  time  of  propagation  is  a  highly  non  linear 
function  of  the  distance  as  illustrated  in  figures  (4,1). 
The  propagation  time  to  any  point  from  the  boundary  here  in 
fact  differs  from  the  change  in  velocity  to  that  point  from 
the  boundary  by  only  a  constant. 

Consider  now  the  case  of  a  plate  being  subducted  into 
the  earth* s  mantle.  The  lower  part  of  the  plate  deforms 
smoothly  and  no  instabilities  occur.  Here  the  lithosphere  is 
a  material  of  lower  melting  point  than  the  surrounding 

mantle  material  (Ringwood  A.  E.  1973)  and  has  a  temperature 
distribution  which  is  fairly  smooth  throughout.  The  initial 
temperature  distribution  is  simply  a  result  of  a  solid 
heated  at  the  lithosphere  asthenosphere  boundary  and  cooled 
at  the  crust  mantle  boundary.  The  temperature  and  stress  at 
the  lower  surface  of  the  plate  then  increase  gradually  as  it 
subducts.  The  subducted  crustal  material  however  is  quite 
different.  It  probably  has  a  lower  conductivity  than  the 
surrounding  mantle  and  is  subducted  with  an  extremely  low 
temperature  relative  to  the  mantle  material.  One  thus  has  a 
large  temperature  difference  at  this  crust  mantle  interface. 
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The  time  for  the  slab  to  reach  thermal  equilibrium  w ith  the 
mantle  when  held  at  a  fixed  position  in  the  mantle  may  be 
estimated  from  t  =p  c  &2/k  to  be  approximately  10a  years. 

Now  consider  the  mimimum  stress  required  to  produce  an 
instability  assuming  constant  stress  throughout  the  slab  and 
an  initial  temperature  equivalent  to  that  of  the  mantle. 
Here  a2=103<*/tM  from  eqn»2.3a.)  where  tM  <10* 3  sec  and  thus 
it  would  take  a  stress  of  10  kilobars  approximately  106 
years  to  produce  an  instability  within  the  crust.  Here  one 
is  also  assuming  the  crust  to  be  thermally  insolated  from 
the  lithosphere.  These  figures  may  not  be  extremely 
reliable*.  but  they  certainly  indicate  that  it  is  not 
reasonable  to  assume  that  instability  can  be  produced  within 
the  slab  by  the  shearing  of  the  mantle.  Here  10  kilobars  is 
far  too  much  shear  for  the  mantle  to  impart  to  the  slab 
since  'the  mantle  motion  induced  by  the  subducting  slab  must 
reduce  its  viscosity  substantially.  If  the  stress  in  the 
slab  is  10  kilobars  however  then  the  stress  in  the  mantle 
must  also  be  10  kilobars.  This  fact  then  yields  a  very 
interesting  result  using  the  same  approximations  as 
previously.  One  now  finds  that  instability  results  in  the 
mantle  in  approximately  102  years  and  within  one  tenth  of  a 
kilometer  of  the  surface  of  the  slab. 

The  previous  analysis  certainly  indicates  that  the 
crust  is  much,  more  stable  than  the  bounding  fluid.  It  thus 
seems  much  more  reasonable  to  consider  the  crust  as 
undergoing  elastic  deformation  and  the  fluid  instability 
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occurring  within  some  boundary  layer  in  the  boundary  fluid. 
This  becomes  an  even  more  attractive  model  when  one 

considers  that  the  crustal  material  is  diffusing  into  the 

mantle  and  thus  inducing  a  compositional  gradient  by  the 
boundary  of  the  slab.  Here  the  lower  melting  point  of  the 
crustal  material  also  means  it  can  support  less  stress  at 
equivalent  temperatures.  If  one  then  also  considers  a  shear 
wave  velocity  gradient  induced  by  a  change  in  chemical 

composition  as  well  as  by  the  change  in  temperature  then 

inertial  motions  due  to  stress  gradients  are  by  far  largest 
at  the  slab  boundary.  The  problem  is  then  whether  the  shear 
heating  induced  by  these  effects  is  large  enough  to  overcome 
the  ability  of  the  material  to  conduct  heat  out  of  this 
region.  If  the  velocity  of  propagation  cf  the  stress  is  many 
orders  of  magnitude  greater  than  the  fluid  velocities  it 
then  appears  reasonable  to  assume  the  stress  is  constant  to 
some  small  depth  6  initially.  Thus  at  6)  =  A  where 
£«a  and  T  ( <$  )  -T  «Tb«  If  the  composition  in  the  material 
being  sheared  was  the  same  throughout  then  the  temperature 
distribution  could  be  written  as 


<K£)  =  (f)(6)  -  2 In 


cosh  ^ (exp (6 ) ) ) 


/  A  (f)(6)  r 

-  /  2  eXP  S—  5J 


Here  the  propagation  time  to  a  distance  E,  from  the  boundary 
is  given  by 
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and  the  velocity  3  may  then  be  written  as 


Now  assume  there  exists  a  depth  <<  <5  such  that 

30«)  *  exp[  (1"^X] 


and 


v  v  (C)  -►  yQ  exp  [  (l-C)xl 


for  E,  <  y 


In  this  region  we  can  observe  how  the  temperat ure  and 
velocity  fields  will  evolve  initially  from  the  class  of 
possible  stationary  solutions  found  in  section  (3.4) « 


The  case  of  interest  in  chapter  3  was  a  temperature 
increase  towards  the  boundary  of  high  viscosity.  In  this 
chapter  the  opposite  case  is  of  intrest,  a  temperature 
decrease  towards  the  boundary  of  high  viscosity.  When  the 
plate  first  subducts  the  heat  loss  from  the  mantle  material 
to  the  upper  plate  boundary  will  be  greater  than  the  effects 
of  shear  heating  in  that  region  because  of  the  very  large 
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temperature  difference  between  the  crust  and  surrounding 
mantle  material.  However  as  the  plate  boundary  approaches 
mantle  temperature  the  effects  of  shear  heating  are  able  to 
dominate.  Under  these  conditions  one  may  obtain  a  maximum 
temperature  in  the  boundary  layer  as  illustrated  in  figures 
(2.5)  and  (2.9).  In  the  presence  of  inertial  motions  caused 
by  the  large  shear  wave  velocity  gradients  in  this  region  a 
sharply  defined  region  with  the  maximum  temperature  must 
develop.  This  is  a  direct  result  of  the  fact  that  inertial 
motion  must  cause  inertial  heating  (see  section  4.2). 

One  is  now  faced  with  a  problem  since  for  deep 
earthquakes  the  failure  is  not  parallel  to  the  subducting 
plate  but  at  a  30°  to  45°  angle  to  plate  subduction.  One 
method  by  which  this  can  be  overcome  is  as  follows.  The 
pressure  gradient  in  the  earth  at  400  km.  depth  is  about 
3.8  x  103  gm/cm 2/see2  (Stacey  F.  D.  1969).  The  temperature 
variations  due  to  shear  heating  cause  large  density 
variations  in  the  boundary  layer.  Here  the  density  of  the 
mantle  material  is  3.8  gm/cm3  (Stacey  F.  D.  1969) ,  while 
4*2  gm/cm3  and  3.5  gm/cm3  are  probably  reasonable  values  for 
the  density  of  the  slab  and  boundary  layer  at  this  depth 
respectively.  This  yields  an  approximately  zero  body  force 
in  the  mantle,  a  force  per  unit  volume  of 
4  x  1  02  gm/cm 2 /sec 2  in  the  slab  directed  towards  the  center 
of  the  earth,  and  a  buoyant  force  per  unit  volume  of 
3  x  1 02  gm/cm2/sec2  in  the  boundary  layer.  It  is  interesting 
to  note  that  these  two  values  are  of  the  same  order  of 
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magnitude,  they  are  however  only  approximations.  The  buoyant 
force  in  the  boundary  layer  for  example  may  vary 
dramatically  with  the  temperature  and  composition  of  the 
material.  Now  assume  that,  the  only  forces  of  importance  in 
the  boundary  layer  are,  the  force  parallel  to  the  plate  due 
to  the  viscous  shearing  of  the  fluid  and  the  buoyant  force 
caused  by  the  heating  of  the  material.  The  axis  of  principal 
shear  at  any  position  in  the  boundary  layer  is  then  oriented 
in  the  direction  of  the  vector  sum  of  these  two  forces. 

This  analysis  may  then  be  summarised  as  follows,  as  the 
boundary  layer  is  sheared  it  is  subjected  to  large  decreases 
in  density  and  is  in  the  neighbourhood  of  melting.  These 
density  changes  cause  compressions!  stresses  and  thus  alter 
the  axis  of  principal  shear.  The  earthquake  must  therefore 
be  a  combination  of  the  relief  of  the  elastic  stress  built 
up  in  the  shearing  of  the  subducting  slab  and  the 
compr essional  stresses  induced  by  temperature  and  pressure 
gradients  in  subduct  ion. 

The  occurence  of  compressions!  earthquakes  thus  seems 
unlikely  at  shallow  depths  for  a  uniformly  moving  plate  due 
to  the  dominance  of  cooling.  It  therefore  seems  that  such  an 
earthquake  would  require  a  large  jerk  to  start  it  on  its  way 
to  fluid  instability.  This  however  is  certainly  not  a 
problem  as  the  plate  motions  are  not  uniform  (Isacks  B.  and 
Molnar  P.  1971,  Spence  W.  1977).  The  frictional  sliding  of 
the  converging  plates  at  the  surface  and  the  tensile 
earthquakes  where  plate  bending  occurs  cause  large  non- 
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uniformities  in  the  motion  of  the  subducting  plates.  These 
compressional  strain  pulses  propagate  down  the  slab  at 
velocities  of  the  order  of  100  km/year  (Spence  W.  1977).  In 
the  following  section  I  will  examine  the  possibility  of 
these  inertial  effects  leading  to  instability. 
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Figure  4c 


The  propagation  time  from  the  boundary  x=1  to 
any  position  in  the  interior  is  given  in 
seconds.  These  times  are  probably  much  smaller 
than  the  propagation  time  for  an  actual  strain- 
rate  pulse  as  the  velocity  of  propagation  is 
assumed  to  be  the  seismic  shear  wave  velocity. 
The  contours  values  for  a  more  slowly 
propagating  plastic  wave  however  should  only 
differ  in  magnitude. 
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1  A  Model  of  Deep,  Earthquakes. 

I  wish  to  consider  here  the  crust  as  an  infinite 
elastic  plate  of  thickness  h  shearing  a  highly  viscous 
fluid.  The  fluid  is  composed  of  crustal  material  at  the 
interface  and  its  composition  changes  to  that  of  mantle 
material  within  a  small  boundary  layer.  The  deformation  of 
the  plate  will  be  assumed  to  be  homogeneous  and  only  in  the 
direction  opposite  to  the  plate  motion.  Since  the  strain  is 
constant  across  the  plate  for  a  homogeneous  deformation  the 
stress  tensor  aik  must  also  be  constant.  If  the  z 
coordinate  is  the  direction  of  plate  motion  and  the  x 
coordinate  is  orthogonal  to  the  plate  surface  then  the  only 
external  force  on  the  plate  is  horizontal  shear  in  the  x 
direction.  For  small  displacements  the  only  non  vanishing 
component  of  the  stress  tensor  is  then  axz. 

\ 

The  stress  for  this  homogeneous  sheared  slab  is  then 
related  to  the  strain  by  (Landau  L.D.  and;  Lifshitz  E.M. 
1959} 

E 

°xz  1+P  exz 

Here  E  is  the  modulus  of  extension  or  Young’s  modulus  and  p, 
the  ratio  of  the  transverse  compression  to  the  longitudinal 
extension,  is  called  Poisson’s  ratio.  The  displacements  are 
thus  related  to  the  stress  by 
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The  velocity  of  slab  subduction  is  of  the  order  10 
cm/year  while  the  velocity  of  plastic  wave  propogation  into 
the  mantle  is  orders  of  magnitude  greater.  Thus  initially 
one  obtains  a  stress  which  is  essentially  constant  to  great 
distances  into  the  surrounding  mantle.  Appendix  C  however 
indicates  that  instability  must  occur  for  any  constant, 
stress  applied  to  a  half  space  and  thus  that  instability 
should  occur  at  a  great  distance  from  the  slab.  This  is  not 
relevant  to  the  problem  at  hand  as  can  be  seen  by  a  simple 
calculation.  If  one  assumes  instability  will  result  and 
simply  assumes  that  it  may  be  expressed  in  terms  of  the 
breakdown  of  the  stationary  stability  criterion  one  then 
obtain  the  following  results.  Instability  at  1000  km  will 
result  for  an  applied  stress  of  approximately  3  bars  and  a 
strain  rate  or  tt  x  10” 16  exp(  $  }  /sec  which  yields  a  plate 
velocity  of  10  cm/year.  The  time  for  this  instability  to 
evolve  assuming  no  heat  conduction  is  then  approximately 
IQio  years  which  is  more  than  twice  the  age  of  the  earth.  In 
any  case  if  one  is  to  accept  the  results  of  section  3.4  that 
50  bars  stress  is  transmitted  by  the  seismic  low  velocity 
zone  then  the  mantle  material  beneath  must  be  capable  of 
imparting  an  even  larger  stress  on  the  subducting  plate. 

By  this  analysis  one  finds  that  for  every  factor  of  10 
increase  in  the  stress  that  the  mantle  is  able  to  offer*  the 
distance  to  instability  decreases  by  a  factor  of  10.  The 
velocity  of  the  plate  in  all  cases  remains  10  cm/year  but 
the  time  to  instability  decreases  in  factors  of  102.  The 
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ability  of  the  mantle  material  to  act  as  a  viscous  fluid  in 
the  neighbourhood  of  the  slab  and  yet  mere  resemble  a  solid f 
beyond  some  depth  Z  in  its  ability  to  resist  flow,  clearly 
results  from  both  geologic  heterogeneity  and  the  fact  that 
the  transition  from  a  highly  viscous  fluid  to  a  solid  is  not 
sharply  defined* 

The  properties  of  these  fluids  may  be  described  as 
follows*  When  the  stress  is  applied  they  undergo  elastic 
deformation.  If  the  deformation  ceases,  shear  stresses 
remain  although  these  are  damped  in  the  course  of  tirae  as 
they  slowly  yield  to  plastic  deformation  (Landau  L.D.  and 
Lifshitz  E.M.  1959)*  The  time  for  stresses  to  be  relieved 
after  the  deformation  stops  is  called  the  “Maxwellian 
relaxation  time*' .  This  relaxation  time  is  much  longer  for 
the  subducting  slab  than  for  the  surrounding  mantle  because 
of  its  higher  viscosity.  It  is  thus  the  mantle  which  first 
yields  to  plastic  deformation  as  the  slab  is  being  subducted 
into  it.  However  the  stress  relief  of  the  mantle  results  in 
shear  heating.  Thus  since  the  fluid  viscosity  decreases  as 
an  exponential  function  of  temperature  the  region  of 
greatest  net  heat  production  behaves  more  and  more  as  fluid 
subjected  to  the  shearing  of  elastic  solids  on  each  side. 
This  results  in  a  region  of  shear  coupling  at  the  crust- 
mantle  interface  similar  to  the  seismic  low  velocity  zone 
beneath  the  plates. 

Now  consider  the  possibility  of  constructing  stable 
solutions  (i.e.  stationary  solutions).  For  a  stationary 
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solution  to  evolve  the  inertial  terms  in  the  stress  strain 
rate  relation  must  vanish  (i.e.  rate  of  stress  relaxation 
must  balance  the  rate  of  stress  application) .  Assume  that 
t-he  subducting  slab  causes  a  zone  of  10  km  thickness  to 
exist  in  a  marginally  stable  state.  This  is  approximately 
the  depth  of  boundary  layer  heating  argued  by  Sydora  (1977) . 
If  we  assume  that  the  normal  mantle  temperature  at  the  depth 
under  consideration  is  200G°K  then  the  stress  in  the  mantle 
is  50  bars.  The  temperature  at  the  boundary  of  this  region 
is  then  2120°C  or  if  we  use  the  stress  strain-rate 
relation  (3«9)  it  becomes  2300°C 

Now  consider  the  boundary  layer  of  composition  change. 
Assume  this  region  has  a  thickness  of  approximately  one 
kilometer.  In  this  region  we  can  approximate  the  heating 
time  by  t  «,=  4x10£S  sec  and  the  cooling  time  by  t  =  4x10J<* 
sec  (section  2.1).  Using  the  analysis  for  a  material  of 
uniform  composition .  Thus  if  the  cooling  is  dominant  over 
the  heating  in  this  approximation  one  obtains  a  temperature 
decrease  at  the  slab  boundary.  But  this  does  not  imply 
stable  solutions  result. 

Instability  may  evolve  in  two  ways.  First  the  crustal 
material  has  a  lower  melting  point  than  the  surrounding 
mantle  material  (Green  T.H.  and  Ringwood  a.e.  1966)  and  one 
obtains  heating  of  the  crust  due  to  thermal  conduction  from 
the  mantle.  The  lower  thermal  conductivity  and  higher 
viscosity  of  the  crust  will  allow  this  boundary  to  evolve 
initially  as  a  boundary  of  thermal  insulation.  However  as 
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the  boundary  heats  it  will  undergo  fluid  deformation  and  due 
to  its  low  melting  point  not  be  able  to  support  the  amount 
of  stress  that  the  mantle  supports.  This  follows  directly 
from  eqn.  (z,4)  and  the  attendant  reduction  in  viscosity 
which  must  occur.  The  fluid  then  must  experience 
accelerations  in  order  to  reduce  the  stress  imparted  to  it. 
If  this  boundary  layer  is  very  thin,  the  onset  of  melt  will 
result  in  larger  accelerations  and  thus  total  melt.  This 
occurs  since  the  applied  stress  remains  greater  even  in  the 
presence  of  partial  melt,  than  this  thin  slab  can  support. 

The  second  way  instability  may  evolve  is  by  obtaining  a 
maximum  temperature  away  from  the  boundary  of  the  subducting 
slab  due  to  shear  heating.  The  viscosity  gradients  in  this 
region  must  result  in  inertial  heating.  Then  as  this 
boundary  layer  is  heated  by  thermal  conduction  one  obtains 
changes  in  the  relaxation  times  of  this  material.  This  in 
turn  causes  stress  gradients  and  thus  accelerations  of  the 
fluid  in  this  region  (sec  4.2}  .  In  this  fashion  one  obtains 
additional  heating  due  to  the  applied  stress  of  the 
subducting  slab  which  is  accentuated  by  the  existence  of 
shear  wave  velocity  gradients.  Thus  once  a  slab  with  a 
concentration  gradient  undergoes  inertial  heating  in  the 
interior  due  to  an  applied  stress  this  in  turn  may  cause 
inertial  motions  which  themselves  become  heat  sources.  This 
means  that  for  a  constant  stress  at  the  boundary,  inertial 
motions  will  occur  in  an  attempt,  to  relieve  excess  heat 
production  at  any  position  in  the  slab.  These  accelerations 
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however  cause  greater  velocity  gradients  at  the  positions 
where  they  occur  and  thus  even  greater  heat  production.  At 
the  slab  boundary  these  accelerations  are  not  able  to  supply 
a  noticeable  reduction  in  stress  due  to  its  high  viscosity 
in  releition  to  the  mantle.  The  problem  then  becomes  whether 
a  real  instability  results  in  which  one  obtains  snapping 
{ stress  relief  due  to  melting  in  a  thin  layer)  or  it  is  able 
to  deform  smoothly  as  the  mantle  material  does. 

The  answer  to  this  problem  depends  on  two  factors.  One 
is  the  amount  of  beating  which  must  take  place  in  order  for 
melting  to  occur  while  the  other  is  how  well  defined  the 
melting  point  is.  If  the  initial  temperature  of  the  material 
is  far  from  melting,  then  difference  in  the  rate  of  heating 
will  cause  a  very  narrow  region  of  melt.  Also  if  the  melting 
point  is  well  defined  then  the  region  which  first  reaches 
this  point  exhibits  a  substantially  lower  viscosity  than  the 
rest  of  the  slab.  It  then  yields  to  almost  the  entire 
deformation  which  is  taking  place  thus  causing  total  melt 
and  a  complete  relief  of  stress. 

Now  consider  the  role  phase  changes  in  the  mantle  and 
and  subducting  slab  may  play  in  this  process.  If  the 
temperature  increases  dramatically  at  the  slab  boundary,  a 
phase  change  will  occur  on  both  sides  of  the  boundary  layer 
prior  to  its  occurrence  within  that  layer.  The  layer  will 
thus  become  narrower  with  depth.  The  material  on  each  side 
will  then  be  able  to  support  more  stress  and  will  have  a 
slower  relaxation  time  than  it  had  before  undergoing  a  phase 
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change.  The  fluid  deformation  carried  by  slab  subduction 
will  thus  be  forced  to  occur  almost  totally  within  the 
boundary  layer.  The  increasing  viscosities  outside  this 
region  will  thus  prevent  a  large  stress  reduction.  Thus  the 
attendant  reduction  in  viscosity  can  easily  cause  the 
stability  criterion  to  be  exceeded  within  the  boundary 
layer.  For  a  viscosity  of  10*  9  poise,  a  thickness  of  1  km 
and  a  stress  of  100  bars  one  obtains  instability  even  for 
the  steady  state  equation  of  motion.  This  however  yields  a 
slab  velocity  of  30  cm/year.  If  the  stress  under  these 
conditions  is  considered  to  increase  to  Ikb  the  stability 
criterion  can  be  exceeded  by  a  viscosity  of  102£  poise  over 
the  1  km  boundary  layer  the  velocity  of  slab  subduction 
necessary  in  this  case  in  the  only  3  cm/year.  This  stress  of 
Ikb,  is  in  fact  approximately  the  stress  obtained  for  slab 
subduction  by  Turcotte  and  Schubert  (1973). 

Finally  consider  the  evidence  that  this  is  in  fact  the 
mechanism  of  deep  earthquakes.  First  the  focal  mechanism  is 
observed  to  be  the  same  as  in  shear  fracture.  In  theory  a 
shear  heating  instability  should  be  indistinguishable  from 
shear  fracture.  Second  this  mechanism  should  cause  volcanic 
activity  following  large  earthquakes.  To  test  correlations 
between  earthquakes  in  the  upper  mantle  and  volcanic 
activities  Blot  ^  1 9 7 7 )  has  attempted  since  1963  to  predict 
volcanic  eruptions  in  the  island  arc  comprising  the  New 
Hebridies,  the  Solomons  and  New  Zealand.  In  the  period  1970 
to  1975  he  forcast  25  eruptions  of  which  20  took  place  with 
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an  accuracy  of  ±15  days*  Other  examples  of  volcanic 
preditions  by  the  method  also  exist  (Grover  J.  1377} .  A 
third  factor  is  that  one  should  expect  an  increase  in  the 
number  of  earthquakes*  with  depth  due  to  greater  mantle 
viscosities  and  with  the  speed  of  subduction.  This  is  also 
observed.  Fourth  the  observed  compressional  earthquakes 
should  define  the  upper  edge  of  the  slab  and  the  orientation 
should  define  the  direction  of  subduction.  These  two 
phenomena  have  also  been  consistently  observed  (Isacks  B. 
and  Holnar  P.  1971)  however  in  the  past  the  earthquakes  have 
been  assumed  to  occur  within  the  upper  surface  of  the  slab. 
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Appendix  A 


An  Analysis  of  the  Variational  Principle  in  Chapter  2 

The  mathematical  approach  presents  here  originates  from 
a  technique  first  used  in  relativistic  particle  dynamics 
(Israel  W .  and  Bailey  I . ,  1975) .  It  is  essentially  an 
attempt  to  insert  in  the  variational  formulation  the  notion 
of  an  irreversible  system.  The  approach  may  be  justified 
rigorously  by  constructing  an  action  functional  simular  to 
the  construction  in  Hamilton*s  principle.  The  fundamental 
difference,  however,  is  that  here  the  action  functional  is 
only  constructed  for  an  infinitesimal  time  step. 
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This  principle  as  applied  to  a  system,  with  a  first 
order  time  derivative  in  the  field  equations,  makes  it  a 
property  of  the  variational  principle  that  the  system  be 
irreversible*  This  is  done  by  defining  the  action  of  the 
system  for  only  an  infinitesimal  time  step  and  constructing 
the  total  action  as  a  sum  of  infinitesimal  actions*  The 
direction  of  the  summation  may  then  be  performed  in  only  the 
future  pointing  time  direction  since  the  initial  conditions 
are  determined  by  the  final  conditions  of  the  previous  time 
step*  It  is  interresting  to  note  that  this  action  functional 
may  also  be  defined  in  terms  of  an  integral  over  a  finite 
time  interval  which  produces  the  equations  of  motion  at 
every  instant.  The  direction  of  evolution  of  the  system  is 
defined  by  the  sign  of  c  which  becomes  the  differential  in 
this  integral.  This  second  type  of  construction  does  not 
lead  itself  to  the  numerical  evaluation  of  the  integral, 
which  is  of  primary  importance  in  the  geophysical 
considerations  for  which  it  is  used  here. 

Tests  on  the  accuracy  of  this  method  have  been  done 
numerically  by  comparing  the  results  obtained  by  Gruntfest 
to  the  solution  of  the  same  problem  with  this  variational 
principle.  The  solutions  obtained  were  identical. 

An  analytical  test  has  also  been  done.  Here  a  semi- 
infinite  slab  initially  at  zero  temperature  has  its  face 
suddenly  raised  to  temperature  TQ .  It  is  assumed  that  T  can 
be  represented  in  the  form:  T=TQ  ( 1 -x/f  ( t)  )  2  where  f  (t)  is 
the  penetration  depth.  Solving  the  variational  principle  for 
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fit)  I  obtain  T=TQ  { l-z/3.  1 6)  2  where  z  =  x/(kt)i/2.  This  is  the 
same  result  as  obtained  by  Djukic  and  Vujanovic  |1971  )  for 
this  problem  and  is  a  good  approximation  to  the  exact 
solution  T  =  T0(1-er£ (z/2)  )  . 
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Appendix  B 


Equations  for  Multiphase  Viscous  Flow  when  Subject  to  a 

Shear  Stress. 

The  equation  of  continuity  for  the  total  mass  of  the 
fluid  is  the  same  as  for  a  one  phase  fluid. 

+  div  (pv)  =0 


Here  the  velocity  is  defined  in  terms  of  the  total  momentum 
of  a  unit  mass  of  fluid.  The  form  of  the  Navier  Stokes 
equation  also  remains  in  the  same  form.  It  however  does 
depend  on  the  concentrations  of  the  different  phases  by 
virtue  of  the  fact  that  viscosity  depends  on  the  relative 
consent rat ions. 


Diffusion  of  the  phases  and  the  solid  body  motion  of 
small  portions  of  the  fluid  both  cause  changes  in  the 
concentrations  of  the  phases  at  different  positions  within 
the  fluid.  The  composition  of  any  given  fluid  element  in  the 
absense  of  diffusion  is  given  by 

dp  _ 


dt  3 1 


~  +  v  .  v  n  =  o 


where  D  is  the  ratio  of  the  mass  of  the  iTH  phase  to  the 
mass  of  the  fluid  in  a  small  volume.  Using  the  equation  of 
continuity  of  total  mass  and  introducing  the  density  of 
diffusion  flux  1  this  equation  becomes 


pc(|—  +  v  «  Vn) 


1  . 

I 


1  18 


The  equation  of  heat  transfer  than  has  an  additional  term 
due  to  diffusion  and  roay  be  written  as 


Here  pn 


pc.  H  _  k  v2d 

a  dt  a 


is  the  chemical 


potential  of  the  i*th 


phase. 


If  the  concentration  of  the  phase  we  are  describing  by 
the  equation  of  continuity  is  small,  as  in  the  case  of  melt 
in  the  ast henosphere,  we  may  express  the  diffusion  flux  as 


i  =  ~pcD  grad  q 

Here  D  the  diffusion  coefficient  and  gives  the  diffusion 


flux  in  this  case 
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Appendix  C 

A  Calculation  of  the  Maximum  Depth  to  which  a 
Stationary  Solution  May  Exist  in  a  Homogeneous  Fluid 
Subjected  to  Constant  Shear. 

Assume  we  have  an  infinite  half  spare.  First  we  wish  to 
consider  the  solutions  for  the  temperature  field  if  the 
stress  is  assumed  constant  for  all  time.  In  this  case  the 
solutions  must  be  stationary. 

2 

k— |  =  -a2  exp  [a  (T-TB)  3 /yQ 
dx" 


a  =  yQ  e  exp  (~4>) 


From  the  temperature  equation  we  obtain 


-  cosh 


-1 


^exp  ((|)m-cJ))  )^j 


*  k 


» 
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exp 
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But  the  boundary  condition  at  x=  00  cannot  be  satisfied 
unless  o  =0.  Wow  assume  we  start  with  a  zero  stress  a  and 
increase  it  infinitely  slowly  so  that  the  stress  can  be 
considered  constant  across  any  thickness  i  .  We  now  wish  to 
find  the  thickness  i  over  which  a  stationary  solution  may 
result.  We  thus  have  the  condition 

$(£)=<$  M~1e2 

The  thickness  i  over  which  a  stationary  solution  exists  is 


In  any  real  physical  situation  however  the  stress  is  not 
increased  infinitely  slowly.  Thus  the  above  relation  given 
an  upper  bound  for  the  length  over  which  a  stable  solution 
may  exist. 


